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1.  Introduction 


The  goal  of  calibrated  spectral  estimation  is  to  determine  sound  pressure 
level  versus  frequency,  at  the  source,  versus  aspect  angle.  Mapping  measured 
spectra  to  the  source  and  knowing  aspect  angle  of  the  source  require  knowledge 
of  the  location  and  heading  of  the  source.  Furthermore,  when  there  are 
multiple  sources,  they  must  be  individually  tracked,  and  their  spectra  must  be 
individually  resolved  to  the  extent  possible.  Consequently,  multi-source 
localization  and  tracking  techniques  comprise  an  important  adjunct  to 
techniques  for  calibrated  spectral  estimation. 

The  localization  objective  is  carried  out  using  acoustic  signals  recorded 
underwater  by  a  hydrophone  array.  In  typical  passive  surveillance  scenarios, 
source  location  is  inferred  from  (1)  the  relative  arrival  times  and  (2)  the 
relative  Doppler  shifts  seen  among  signals  received  at  the  various  sensors 
[26]. 


The  relative  arrival  time  between  two  sensors  is  called  the  time- 
difference  of  arrival  (TDOA),  and  a  large  literature  on  delay  estimation  has 
developed  toward  optimal  and  time-varying  estimation  of  TDOAs  (for  example,  a 
special  issue  on  this  topic  appears  in  [21]).  Delay  estimation  can  also  be 
applied  to  the  echoes  seen  at  a  single  sensor  due  to  multipath;  in  effect, 
multipath  reception  can  provide  an  extra  "virtual  sensor"  located  at  the 
mirror  image  of  a  real  sensor  in  the  reflecting  surface. 

Estimation  of  differential  Doppler  shift  is  not  as  well  studied  as  the 
delay  estimation  problem,  and  so  it  will  probably  remain  a  fertile  research 
area  for  some  time.  Unlike  TDOA  estimation  which  requires  a  broadband  signal 
for  its  accurate  measurement,  differential  Doppler  can  be  measured  in  either 
the  narrowband  or  broadband  case.  The  narrowband  case  is  very 
straightforward,  involving  only  ratios  of  FFT  peak  frequencies,  while 
broadband  spectral  cross-correlation  appears  to  offer  significant 
opportunities  for  improvement.  There  seems  to  be  no  mature  treatment  of 
optimal  broadband  Doppler  correlation  such  as  exists  for  delay-based 
correlation.  An  advantage  of  localization  based  on  Doppler  is  that  Doppler 
shifts  can  be  reliably  measured  at  much  greater  range  than  TDDA  under  typical 


ci rcumstances. 


When  there  is  more  than  one  source,  localization  is  typically 
accomplished  by  repeated  application  of  the  single-source  methods.  In  the 
case  of  measuring  TDOAs,  each  source  produces  its  own  secondary  peak  in  the 
cross-correlation  function  between  two  sensors;  if  the  peaks  are  resolved 
(i.e.,  the  signal  bandwidths  are  sufficiently  large  and  the  TOOA  times  are 
sufficiently  separated),  then  there  is  no  problem  obtaining  N  TDOAs  per  sensor 
in  the  case  of  N  sources.  However,  it  is  then  necessary  to  try,  in  principle, 
all  possible  associations  of  TDOA  to  source;  that  is,  N  groups  must  be  chosen, 
each  consisting  of  one  TDOA  from  each  sensor  pair,  and  each  group  must  have 
TDOAs  arising  from  one  source  only.  While  not  imposing  any  fundamental 
barriers,  the  association  problem  can  vastly  increase  the  computational  burden 
in  the  multi-source  case  relative  to  the  single-source  case. 

1.1  The  SPICE  Project 


The  purpose  of  the  SPICE  project  at  SCT  has  been  to  develop  new 
techniques  which  aid  in  the  detection  and  localization  of  multiple  underwater 
sources.  In  this  section,  the  major  components  of  the  project  are  summarized. 

The  first  phase  of  the  project  was  devoted  to  developing  new  methods  for 
high-resolution  spectrum  analysis.  Improvements  in  spectral  resolution  (given 
fixed  statistical  stability)  translate  directly  to  greater  tracking 
accuracy.  In  many  ways,  this  phase  of  the  project  was  a  follow-up  effort  to  a 
prior  project  entitled  "Multi -Target  Tracking  Studies"  (MTS)  which  started  in 
1979  (cf.  SCT  Reports  5334-01,02).  A  mature  treatment  of  the  high-resolution 
techniques  is  presented  in  SCT  Report  5498-04.  See  also  SCT  Report  5465-02. 
All  SCT  Reports  mentioned  in  this  report  should  be  requested  from  the  Adaptive 
Systems  Department  of  the  Advanced  Technology  Division  of  SCT. 


The  second  phase  of  the  SPICE  project  was  devoted  to  developing, 
simulating,  and  evaluating  a  novel  multisource  tracking  algorithm  devised  by 
Dr.  Benjamin  Friedlander.  For  source  detection  in  a  plane,  the  technique  j$[ 

involves  finding  intersections  (one  per  source)  of  five  or  more  hyperplanes  in 
5-space.  The  multisource  tracking  algorithm  is  fully  described  in  SCT  Report 

V 
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The  third  and  final  phase  of  the  SPICE  project  was  to  develop  more 
traditional  sonar  signal  processing  facilities  based  on  nonparametric  spectral 
estimation  using  an  FPS  array  processor.  Estimates  of  intersensor  time  delay, 
multipath  time  delay,  intersensor  differential  Doppler,  and  absolute  Doppler 
are  all  derived  in  one  way  or  another  from  the  FFT-based  spectrum  estimate. 
Both  narrowband  and  broadband  Doppler  estimation  techniques  have  been 
developed.  Differential  and  absolute  Doppler  estimates  have  been  produced 
from  live  sonar  recordings.  The  high-resolution  spectrum  analysis  program  and 
the  FFT-based  power  spectrum  estimator  are  documented  in  SCT  Report  5466- 
06B.  The  differential  delay/Doppler  techniques  are  described  only  in  this 
document. 

Side  investigations  carried  out  on  the  SPICE  project  included  a  new 
technique  for  delay  estimation  (SCT  Report  5466-03)  and  the  application  of 
this  method  to  the  tracking  of  multipath  delay  (SCT  Report  5466-04).  Also,  a 
constrained  adaptive  notch  filter  was  developed  to  enable  the  real-time 
elimination  of  unwanted  sinusoidal  components  ("Analysis  and  Performance 
Evaluation  of  an  Adaptive  Notch  Filter"  by  B.  Friedlander  and  J.  0.  Smith, 
SCT  ATD/ASD,  1982). 

The  present  final  report  for  the  SPICE  project  includes  all  material  not 
appearing  in  prior  reports,  and  provides  a  list  of  all  reports  and 
publications  deriving  from  the  SPICE  project. 

2.  Guide  to  the  Appendices 

The  appendices  included  here  discuss  the  localization  and  tracking  of 
sources  using  TOGA  and  differential  doppler  measurements.  Improvements  to 
currently  existing  systems  are  outlined,  and  several  new  computationally 
efficient  localization  and  tracking  methods  are  presented.  The  associated 
problem  of  measuring  and  tracking  TDOA  and  differential  Doppler  values  is  also 
covered.  Finally,  an  appendix  discussing  an  interactive  signal  processing 
development  environment  is  included.  Below,  brief  descriptions  of  the 
appendices  are  given. 


Appendix  A  discusses  ideas  for  extending  the  Omni -Tracking  System  (OTS) 
to  handle  multiple  sources  more  effectively.  A  wide  range  of  alternatives  is 
explored,  from  ways  to  use  the  existing  system  unmodified  to  a  complete 
replacement  of  OTS  by  a  super-powerful ,  model -based,  track  finder. 

Appendix  B  describes  a  fast  track  solver  devised  by  Dr.  Friedlander.  It 
uses  simple  features  of  "S-curves"  measured  in  an  absolute-Ooppler  path  to 
compute  track  parameters.  While  measurement  of  the  required  features  (such  as 
minimum/maximum  long-range  Doppler)  depends  on  a  high  signal -to-noise  ratio, 
the  track  parameters  can  be  found  very  quickly  relative  to  more  commonly  used 
methods. 

Appendix  C  documents  various  signal  measurements  which  are  useful  for 
multisource  localization.  Most  of  these  are  currently  in  use,  and  the 
appendix  serves  as  a  review  of  some  basic  signal  processing  fundamentals.  In 
addition,  the  autocorrel ogram  and  its  use  in  localization  (described  in 
Appendices  A  and  B)  are  believed  to  be  new  ideas. 

Appendix  D  describes  newly  developed  Spherical  Interpolation  method  for 
localizing  a  single  source  based  on  time-di fference-of-arri val  (TDOA) 
measurements.  The  method  is  very  fast  computationally  and  surprisingly 
accurate,  coming  close  to  the  Cramer-Rao  lower  bound  for  unbiased  estimators, 
as  shown  in  Appendix  E. 

Appendix  F  discusses  track  parameter  estimation  from  multipath 
information.  Mathematical  techniques  analogous  to  those  used  to  develop  the 
Spherical  Interpolation  method  are  used. 

Appendix  G  describes  the  kind  of  interactive  development  environment  we 
feel  is  important  to  provide  the  underwater  surveillance  engineer.  Its  major 
feature  is  malleability  allowing  all  signal  processing  tools  to  be  used  in 
interactive  exploratory  analysis  with  a  complete  set  of  signal  display 
features  and  display  controls.  The  reason  for  choosing  an  interactive 
programming  environment  is  that  one  cannot  know  in  advance  what  signal 
processing  procedures  are  going  to  be  most  effective  for  localizing  a  given 
source.  Instead,  the  measured  data  must  be  analyzed  interactively  for 


features  which  will  give  the  source  away  to  an  automated  analysis.  Another 
important  feature  of  the  proposed  system  is  the  ability  to  record  interactive 
work  for  later  compilation  into  an  automated  procedure.  We  believe  such  a 
system  can  increase  the  effectiveness  per  man-hour  of  localization  software 
development  by  a  large  factor. 

3.  Conclusions 

Several  avenues  for  increasing  the  effectiveness  of  multisource 
localization  have  been  explored  in  this  project.  Here  we  will  indicate  what 
appear  to  be  the  most  promising  directions  for  future  development,  based  on 
our  results. 

We  feel  that  the  most  powerful  approach  to  multisource  localization  is 
the  dynamic  programming  system  which  explicitly  evaluates  the  likelihood  of  a 
source  track  by  comparing  actual  sensor  measurements  to  synthetic  measurements 
generated  by  a  model  of  the  source  moving  along  the  hypothesized  track.  Such 
a  framework  provides  for  incorporation  of  all  a  priori  knowledge  about  the 
source  and  the  scenario.  This  method  is  discussed  briefly  in  Appendix  A  and 
to  a  greater  extent  in  Appendix  F. 

As  discussed  previously,  relative  time  delay  and  relative  Doppler  shift 
between  sensors  provide  the  basis  for  present  localization  systems.  As  stated 
earlier,  time-delay  estimation  is  a  relatively  mature  field,  while  Doppler 
estimation,  especially  in  the  broadband  case,  has  not  been  fully  analyzed  in 
the  literature.  Specifically,  broadband  Doppler  correlation  methods  (as 
discussed  in  Appendices  C  and  G)  need  to  be  analyzed  in  a  manner  analogous  to 
what  Van  Trees  [25]  has  done  for  more  conventional  correlation  methods. 
Additionally,  the  problem  of  converting  instantaneous  TDOA  and  Doppler 
measurements  into  source  location  and  velocity  estimates  has  not  been 
thoroughly  studied.  In  particular,  currently,  there  appears  to  be  only  a  few 
computationally  efficient,  accurate  source  location  estimators  from  TDOA 
measurements  (see  appendices  D  and  E),  and  no  computationally  efficient, 
accurate  estimators  of  source  location  and  velocity  from  instantaneous  doppler 
or  differential  doppler  measurements,  although  there  are  approximate  solutions 
to  the  2-dimensional  special  case  of  a  source  far  away  from  a  linear  array 


£ 
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The  very  effective  Spherical  Interpolation  method  for  localizing  a  source  „ 
from  TDOA  measurements  (Appendices  D,  E,  and  F)  warrants  further  study.  In 
particular,  it  seems  quite  possible  that  the  same  basic  approach  can  be 

r  • 

applied  to  localization  from  differential  Doppler  measurements.  If  this  is  2 

true,  and  if  the  Doppler  case  is  as  close  to  optimal  as  is  the  TDOA  case, 
current  practice  in  passive  underwater  surveillance  will  probably  undergo  a 

y 

major  revision. 
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Appendix  A  -  Extending  the  OTS  System  to  Multiple  Sources 


The  Omni -Tracking  System  (OTS)  [1]  carries  out  the  localization  objective 
for  a  single  source.  The  OTS  system  estimates  the  track  of  a  single 
maneuvering  source  based  on  Doppler  cross-correlograms,  delay  cross- 
correlograms,  or  both.  (The  delay  correlogram  yields  a  time-difference  of 
arrival  (TDOA)  versus  time  between  each  pair  of  sensors.) 

This  appendix  discusses  alternatives  for  extending  the  operation  of  OTS 
to  multiple  sources.  Familiarity  with  the  OTS  is  assumed. 

A.l  Current  Use  of  the  OTS  System  in  the  Presence  of  Multiple  Sources 


When  there  are  multiple  sources,  the  operator  must  correctly  associate 
correlogram  peaks  across  time  to  obtain  the  differential  Doppler  or  delay 
corresponding  to  the  same  single  source  in  each  sensor  pairing.  Such  an 
association  can  be  difficult,  if  not  impossible,  for  an  operator  using  only  a 
display  of  correlogram  peaks  versus  time. 

BEARTRK  is  the  track-produci ng  component  of  OTS.  When  the  operator  passes 
the  delay/Doppler  information  for  a  single  source  to  BEARTRK,  a  nonlinear 
optimization  is  performed  with  respect  to  the  track  parameters.  Because  the 
optimization  is  nonlinear,  convergence  to  the  best  estimate  depends  on 
initial ization  sufficiently  close  to  the  true  track  parameters.  The  use  of  a 
gradient/Newton  descent  method  indicates  that  the  problem  is  too  large  for 
exhaustive  search.  All  such  nonlinear  optimization  procedures  carry  the  risk 
of  convergence  to  a  false  local  minimum. 

The  ADEC  program  [2]  is  used  to  prepare  line  tracks  through  the 
individual  sensor  spectrograms.  An  operator  manually  prunes  the  lines, 
removing  unwanted  lines,  filling  in  partially  missing  lines,  and  labeling  the 
image  of  each  source  line  in  each  receiver.  The  frequency-versus-time  tracks 
so  obtained  are  then  divided  into  all  possible  pairings  (point  by  point)  to 
produce  differential  Doppler  estimates.  To  increase  robustness,  the  Doppler 
ratios  for  several  lines  may  be  averaged  together  (since  they  are  assumed  to 
emanate  from  a  single  source).  The  Doppler  ratios  are  then  fed  to  BEARTRK  for 


localization  in  the  usual  way. 


ADEC  is  a  non-coherent  tracker  for  sinusoidal  lines  in  noise,  and  it 
performs  line  tracking  in  a  manner  somewhat  similar  to  the  coherent  tracker 
MAPLE  [3,4].  ADEC  operates  by  tracking  peaks  through  time  in  the  normalized 
spectrogram,  or  LOFARGRAM.  A  running  estimate  of  frequency  f^(t)  and 
frequency  rate  d^(t)  is  maintained  for  each  line  being  tracked.  The  smoothed 
frequency  and  rate  estimates  are  used  to  predict  the  peak  frequency  for  the 
next  time  step,  and  this  prediction  defines  the  center-frequency  of  a  search 
window  within  which  the  peak  of  the  power  spectrum  is  found.  The  window  width 
is  adaptive  based  on  amplitude.  The  peak  amplitudes  are  used  to  compute  a 
running  figure  of  merit  for  the  tracked  line  related  to  the  likelihood 
function  for  a  sinusoid  in  white  noise.  The  behavior  of  the  incremental 
likelihood  variable  with  respect  to  various  thresholds  determines  when  a  line 
is  detected,  tracked,  and  declared  finished. 

The  ADEC  output  is  generated  from  the  spectrogram  at  every  sensor.  For 
example,  the  lines  in  sensor  2  could  be  labeled  as  A2(t),  B2(t),  C2(t),  and  so 
on.  The  operator  then  decides  which  lines  in  each  spectrogram  are  from  the 
same  original  line,  and  Doppler  ratios  of  the  form  A2(t)/Al(t)  are  formed. 

The  Doppler  ratio  tracks  are  matched  in  BEARTRK  by  searching  for  the  track- 
parameters  which  best  predict  the  observed  Doppler  ratios. 

A. 2  Block-Exhaustive  Search  in  BEARTRK 

A  straightforward  extension  of  the  OTS  system  to  multiple  sources  would 
be  to  accept  all  of  the  correlogram  peaks  produced  by  CORAN  (with  little  or  no 
pruning  by  the  operator)  and  to  fit  multiple  tracks  to  these  peaks.  If 
BEARTRK  would  exhaustively  search  the  parameter  space,  this  extension  would  be 
straightforward.  However,  because  the  track  estimation  is  based  on  gradient 
descent,  it  is  typically  not  possible  for  the  algorithm  to  explore  different 
associations  of  correlogram  peaks  to  sources.  (To  do  so  would  normally 
require  moving  uphill  on  the  error  surface.)  A  mulfisource  version  of  BEARTRK 
could  be  written  which  would  combine  coarse  exhaustive  search  followed  by 
iterative  fine  tuning  by  gradient  descent.  The  result  would  be  several 
locally  optimal  tracks  among  which  t'ne  true  tracks  are  assumed  to  lie.  Only 
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BEARTRK  need  be  modified  using  this  approach  to  OTS  extension. 


A. 3  Exhaustive  Search  over  CQRAN-Peak  Associations  (OASY/PREP). 

Instead  of  developing  an  exhaustive-search  version  of  BEARTRK,  the 
association  of  CORAN  peaks  can  be  automated.  In  other  words,  the  normally 
manual  association  of  correlogram  peaks  across  time  (to  arrive  at  a  path 
corresponding  to  one  source)  can  be  carried  out  automatically.  To  try  various 
associations  of  correlogram  peaks  in  the  multisource  case,  BEARTRK,  as  is, 
could  be  run  in  several  passes,  each  of  which  is  initialized  to  a  different 
sequence  of  correlogram  peaks  versus  time.  Each  pass  of  BEARTRK  would  produce 
a  track  estimate  and  a  measure  of  fit.  After  all  passes  are  complete,  the 
tracks  can  be  sorted  by  measure  of  fit,  and  the  best  ones  chosen  as  track 
estimates. 

An  operator  can  greatly  reduce  the  number  of  CORAN  peak  association 
tracks  to  be  tried.  In  the  automated  case,  some  search  heuristic  reduction 
analogous  to  that  performed  by  the  operator  is  desirable.  As  an  example,  if 
there  are  5  correlogram  peaks  per  time  step  and  30  time  steps,  BEARTRK  would 

on  ?i 

be  called  on  the  order  of  5  *  10  times  if  all  possible  associations  were 

to  be  tried.  If  BEARTRK  could  produce  a  track  estimate  in  one  second,  this 
process  would  take  3Q  trillion  years.  Thus,  not  only  is  exhaustive  search  of 
the  parameter  space  prohibitive  for  individual  tracks,  it  is  also  beyond 
present  computing  technology  to  try  all  possible  associations  of  peaks  to 
correlogram  lines.  Clearly,  it  is  necessary  to  prune  the  associations  and 
pass  only  a  few  "reasonable’1  possibilities  to  BEARTRK.  One  method  for  cutting 
down  the  search  possibilities  is  to  apply  the  Viterbi  Algorithm  [5]. 

A. 4  The  Block-Viterbi  Algorithm  (BVA). 

The  standard  method  for  finding  multiple  paths  through  time  in  a 
computationally  feasible  way  is  called  dynamic  programming  or  the  Viterbi 
algorithm  [5].  The  idea,  in  the  first-order  case,  is  to  evaluate  all  paths  in 
parallel  through  time,  saving  only  the  best  paths  at  each  time  step  for 
extension  to  the  next.  In  the  5-peak,  30-step  example  above,  there  would  be 
on  the  order  of  5  •  30  =  150  path  "evaluations"  per  source  in  contrast  to 


lO^l.  (A  path  is  now  defined  as  a  connection  of  correlogram  peaks  through 
time  for  a  given  sensor  pair  -  this  is  not  to  be  confused  with  a  track  which 
still  refers  to  the  position  and  velocity  history  of  a  source  through  time.) 

Unfortunately,  the  first-order  Viterbi  algorithm  is  not  immediately 
applicable  to  the  line-association  problem.  To  see  this,  consider  that  the 
first  time  BEARTRK  (or  suitable  replacement)  is  called  to  fit  a  track  to  the 
first  correlogram  peak,  it  has  only  one  point  to  work  with  from  each  sensor 
pair.  The  track  may  be  underdetermined.  Therefore,  we  must  be  able  to  go  to 
a  higher  order  Viterbi  algorithm  which  looks  at  several  time  steps  before 
making  a  decision  on  path  extension.  A  Kth-order  Viterbi  algorithm 
exhaustively  evaluates  K  time  steps  per  path  extension  (placing  the  number  of 
path  evaluations  somewhere  between  150  and  lO^1).  That  is,  the  first  K  time 
steps  are  tested  exhaustively,  the  winning  Ns  paths  are  extended  by  K  points 
in  the  same  way,  and  so  on.  (Ns  is  the  number  of  paths  "kept  alive"  at  each 
step.)  In  the  above  example,  if  the  order  is  K=3  and  Ns=l  (one  source),  we 
get  53(30/3)=1250  path  evaluations  for  one  source. 

More  generally,  for  Np  peaks  per  correlogram,  Nt  time  steps,  Ns  sources, 
and  order  K  time  steps  per  path  extension,  the  "bl ock-Viterbi "  algorithm  (BVA) 
requires  at  most 
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Typically,  this  number  refers  to  the  number  of  incremental  path  evaluations.  v 

The  cost  of  each  "incremental  evaluation"  is  the  cost  of  updating  likelihood 
of  the  whole  path  to  include  K  new  steps  forward.  By  using  log  likelihood  to  ^ 
evaluate  the  path,  these  updates  are  additive. 

In  practice,  the  number  Nk  can  be  replaced  by  Nk  where  N  is  the  hi 

p  mm 

maximum  number  of  correlation  lags  a  path  can  change  in  one  time  step 

(typically  only  one  or  two).  -S 


There  is  still  the  problem  that,  far  from  CPA,  the  track  estimate  is  very  .V 

Li 

noisy.  This  suggests  starting  the  BVA  at  the  average  CPA  time  and  working 

outward  into  the  past  and  future.  This  can  be  an  important  refinement  because  r-. 
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the  Viterbi  algorithm  does  not  reconsider  its  early  decisions. 

A. 5  Lag-Limited  Exhaustive  Association 

Since  a  realistic  correlation  path  cannot  change  very  quickly  from  time 
step  to  time  step  (say  less  than  Nm  correlation  lags  per  time  step), 
exhaustive  search  by  dynamic  programming  can  be  made  much  more  efficient  when 
this  constraint  is  employed. 

Wolcin  [4]  proposes  an  effective  strategy  for  incorporating  a  time-rate- 
of-change  constraint  in  the  context  of  using  MAPLE  [3,4]  to  track  slowly 
changing  sinusoidal  frequencies  in  a  spectrum.  The  idea  is  based  on  the 
observation  that  each  possible  initial  path  point  expands  into  a  small  "wedge" 
of  possible  paths.  For  example,  if  N  =  1  (defined  above),  then  the  number 
of  points  which  must  be  considered  at  each  time  step  grows  as  1,  3,  5,  7,  and 
so  on.  At  each  time  step,  only  two  more  points  are  within  reach  of  a  path 
from  a  given  point  at  time  0. 

Wolcin's  idea  is  to  choose  a  set  of  starting  points  which  are  close 
enough  together  so  that  their  "possible-path  wedges"  intersect  half  way  across 
the  total  time  span.  For  example,  if  there  were  30  time  steps  and  =  1  as 
before,  then  path  starting  points  would  be  chosen  every  30  correlation  lags, 
because  the  two  wedges  emanating  from  starting  points  30  lags  apart  will 
overlap  each  other  after  15  time  steps.  The  supposition  is  that  a  true  path 
will  "capture"  one  of  the  expanding  wedges  of  exhaustive  search;  therefore,  it 
is  not  necessary  to  try  all  possible  starting  points  to  detect  the  presence 
(and  latter  trajectory)  of  a  path.  However,  to  obtain  the  initial  trajectory 
of  the  path,  the  Viterbi  algorithm  is  run  backwards  in  time  through  the 
correlogram  beginning  at  the  final  point  of  the  captured  path.  As  a  final 
refinement,  Wolcin  recommends  pursuing  the  path  forward  once  again  to  verify 
that  the  best  path  from  the  newly  found  correct  starting  point  in  fact  ends  on 
the  point  used  for  the  backward  search.  If  the  third  pass  diverges  from  the 
latter  path  segment  picked  up  on  the  first  pass,  the  path  is  rejected. 

Detected  paths  are  then  removed  by  zeroing  the  correlation  bins  through 
which  the  path  moved.  This  zeroing  must  be  done  in  a  manner  matched  to  the 


analysis  used  in  obtaining  the  correlogram.  If  the  signal  or  analysis 
bandwidth  is  small,  for  example,  then  the  correlogram  paths  are  characterized 
by  wide  peaks  (covering  multiple  correlation  lags).  In  this  situation,  a 
simple  zeroing  of  correlogram  along  the  path  lags  is  not  efficient  because 
large  sharp  ridges  will  be  left  on  either  side  of  the  zeroed  path  which  can  be 
interpreted  as  two  more  paths  running  in  parallel  with  the  first.  In  general, 
the  path  zeroing  should  consist  of  multiplication  by  an  appropriate  "window" 
function  whose  shape  equals  the  inverse  of  the  expected  cross-correlation  peak 
shape  between  signals  received  at  two  different  sensors. 


A. 6  Path  Association  before  Track  Computation 

All  multisource  tracking  alternatives  require  a  track-fitting  procedure 
(such  as  8EARTRK)  to  be  called  for  each  path  evaluation.  A  potential  speed-up 
is  to  first  determine  a  set  of  all  feasible  paths  through  the  correlogram 
peaks.  The  set  of  feasible  paths  is  all  possible  paths  minus  those  deemed 
physically  unrealistic.  This  automatic  pruning  step  can  greatly  reduce  the 
size  of  the  search  space.  Note  that  feasibility  constraints  are  easily 
incorporated  into  the  81 ock-Viterbi  Algorithm.  The  lag-limited  search 
described  in  the  previous  subsection  is  one  special  case  of  feasibility 
constraints. 

A. 7  Dynamic  Programming  through  Unprocessed  Correlograms. 

Instead  of  working  with  the  output  of  CORAN,  which  consists  of  up  to  five 
correlogram  peaks  for  each  time  step,  a  multiple-path-finding  algorithm  can  be 
applied  directly  to  the  delay  or  Doppler  correlogram  for  each  sensor  pair. 

This  allows  peaks  to  be  tracked  across  time  frames  in  which  they  are 
momentarily  below  the  five  largest  peaks  of  the  correlogram.  In  other  words, 
once  a  path  is  under  way,  the  nearest  local  maximum  (of  sufficient  amplitude) 
in  the  correlogram  can  be  taken  as  the  continuation  of  the  path,  rather  than 
depending  on  the  path  continuation  lying  among  the  five  largest  correlogram 
peaks.  If  there  are  time  frames  in  practice  where  a  correlogram  peak  is 
missing,  direct  peak  tracking  in  the  correlogram  should  be  considered. 
Similarly,  dynamic  programming  can  be  applied  to  the  spectrogram  (as  in  ADEC) 
or  auto-Doppler  correlogram  (cf.  Appendix  C)  for  the  purpose  of  measuring 


absolute  frequency  or  Doppler  versus  time. 


Sinusoidal  signal  components  can  be  tracked  by  a  coherent  estimator,  as 
is  done  by  MAPLE  [3,4].  Coherent  processing  yields  a  3dB  improvement  in  the 
signal -to-noise  ratio  (SNR),  relative  to  non-coherent  processing,  for  every 
doubling  of  the  integration  time. 

The  frequency  estimate  produced  by  a  non-coherent  frequency  estimator  is 
the  (interpolated)  location  of  the  peak  of  the  power-spectral  density  (PSD) 
estimate.  The  PSD  estimate  is  normally  computed  as  a  time-average  of  the 
magnitude-squared  FFT's  of  successive  time  frames  [6],  The  FFT  length 
determines  the  height  of  a  sinusoidal  peak  above  the  noise  floor. 

A  coherent  frequency  estimator  conceptually  maximizes  the  magnitude  of 
the  inner  product  between  the  entire  time  signal  and  a  sinusoid  oscillating  at 
the  estimated  frequency.  The  effective  height  of  the  sinusoidal  component 
above  the  noise  floor  grows  with  the  length  of  signal  processed  (3dB  for  every 
doubling  of  observation  time).  For  the  case  of  a  single  sinusoid  in  white 
Gaussian  noise,  this  estimate  coincides  with  the  maximum-likelihood  estimate 
(also  the  minimum-variance  estimate  in  the  non-Gaussian  white  noise  case)  [7]. 
The  use  of  a  coherent  estimate  requires  enough  data  to  reach  an  effective 
si gnal -to-noi se  ratio  per  bin  much  greater  than  zero.  (The  unaveraged 
magnitude-squared  of  the  Fourier  transform  of  Gaussian  white  noise  has  a  (chi- 
squared)  standard  deviation  equal  to  its  mean  [fi].) 

In  the  case  of  a  time-varying  frequency  (due  to  Doppler),  the  maximum- 
likelihood  estimate  requires  searching  over  all  possible  time-varying 
frequency  histories  and  again  maximizing  the  coherent  inner  product  with  the 
entire  time  signal.  The  MAPLE  algorithm  provides  an  approximate  solution  to 
this  problem  by  using  dynamic  programming  to  incrementally  maximize  the  inner 
product  of  the  time  signal  with  a  sinusoid  having  piecewise  linear  frequency 
variation.  In  addition,  MAPLE  allows  a  penalty  to  be  placed  on  frequency 
change,  and  it  can  be  confined  to  search  paths  only  in  a  small  frequency 
interval  (e.g.,  the  known  sinusoidal  frequency  plus  and  minus  the  maximum 


expected  noppler  shift).  It  is  straightforward  to  include  all  realistic 
physical  constraints  associated  with  a  Doppler-shifted  sinusoidal  underwater 
source.  The  highly  restricted  variation  in  this  case  makes  MAPLE  not  nearly 
as  computationally  expensive  as  it  is  when  searching  a  wide  frequency  region 
and/or  fast  frequency  variation.  The  MAPLE  algorithm  in  its  present  form  can 
be  used  to  greatly  extend  the  range  over  which  the  Doppler-shift  of  individual 
sinusoidal  lines  can  be  tracked. 

An  important  feature,  not  employed  by  MAPLE  and  normally  employed  by 
Doppler  processing  within  OTS,  is  the  total  spectral  correlation  used  in  the 
Doppler  correlogram.  The  Doppler  correlogram  is  non-coherent ,  but  uses  the 
whole  spectrum  in  its  Doppler  estimate;  MAPLE  is  coherent,  but  looks  only  at 
one  frequency  at  a  time.  (In  cases  of  high  SNR  and  multiple  sources,  MAPLE 
has  the  advantage  of  not  producing  "cross-terms"  associated  with  spurious 
secondary  correlations  between  the  spectra  of  two  different  sources.)  To 
recover  the  whole-spectrum  advantage  of  the  Doppler  correlogram  in  a  coherent 
tracker  such  as  MAPLE,  an  extension  is  necessary:  MAPLE  can  be  modified  to 
simultaneously  track  a  set  of  spectral  lines,  all  derived  from  a  single 
Doppler  track.  Such  a  version  of  MAPLE  would  be  exceedingly  effective  in 
cases  where  the  true  sinusoidal  frequencies  are  known.  In  operation,  a  single 
value  of  Doppler  (suitably  constrained  in  its  time  variation  and  extent)  would 
be  optimized  using  dynamic  programming  so  as  to  maximize  the  inner  product  of 
the  time  signal  with  the  sum  of  Doppler-shifted  sinusoids. 

A. 9  Dynamic  Programming  through  the  Track  Parameter  Space 


The  ideal  application  of  dynamic  programming  is  to  directly  evaluate  the 
likelihood  of  the  sensor  measurements  as  a  function  of  time  through  all 
possible  source  tracks.  Searching  directly  through  possible  source  tracks  as 
opposed  to  correlogram  peak  tracks  or  spectrogram  peak  tracks  allows  maximum 
use  of  a  priori  knowledge  regarding  the  physical  constraints  and 
characteristics.  Also,  the  likelihood  computation  is  being  applied  to  the 
actual  quantity  of  interest  -  the  track  -  rather  than  some  indirect 
manifestation  of  the  track  such  as  its  correlogram/spectrogram  image  in  sensor 
pairs.  Finally,  when  the  track  space  is  searched  directly,  the  track 
parameter  estimation  function  of  BEARTRK  is  fully  absorbed. 
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The  difference  between  this  approach  and  that  discussed  in  subsection  A. 2 
is  in  the  definition  of  sensor  measurements  which  the  track  estimate  must 
predict.  Use  of  the  raw  sensor  signals  (or  successive  complex  spectra  as  used 
by  MAPLE)  instead  of  CORAN  peaks  allows  incorporation  of  all  known 
envi ronmental  data  such  as  signal -to-noi se  ratio,  noise-floor  shape, 
interference  activity,  and  measurement  stability. 

Full-scale  dynamic  programming  techniques  for  multisource  tracking  are 
described  further  in  section  5  of  Appendix  G. 

A. 10  The  Auto-Ooppler  Correlogram 

The  OTS  system  bases  its  track  estimates  on  differential  Doppler  and/or 
delay  between  sensors.  It  seems  worthwhile  to  consider  other  sources  of 
information  which  may  provide  better  track  estimates  even  in  the  single-source 
case.  In  Appendix  C,  the  Auto-Ooppler  correlogram  is  described.  Essentially, 
instead  of  cross-correlating  the  spectrum  in  one  sensor  with  that  in  another, 
the  short-time  PSD  estimate  in  a  single  sensor  is  cross-correlated  against  the 
PSD  in  the  same  sensor  at  a  fixed  time. 

We  have  found  that  for  individual  spectral  lines  (visible  in  the 
spectrogram  of  PSD  versus  time)  the  time  of  CPA  can  be  estimated  as  the  time 
of  maximum  intensity  for  the  line.  By  choosing  the  CPA  time  of  a  particular 
line  as  the  reference  time  for  the  auto-Doppler  correlogram,  the  auto-Doppler 
correlogram  becomes  largely  "matched”  to  the  spectrum  of  the  source  emitting 
the  reference  line,  and  the  peak  track  in  the  auto-Doppler  correlogram  is 
likely  to  provide  a  good  estimate  of  the  absolute  Doppler  shift  versus  time 
for  the  source  whose  CPA  was  determined.  This  procedure  can  be  repeated  for 
each  distinct  time  of  CPA  for  each  line  in  each  sensor.  Presumably,  these 
auto-Doppler  correlograms  would  produce  a  set  of  absolute  Doppler  estimates 
for  each  source  present.  This  information  can  be  combined  with  differential 
Ooppler  estimates  to  sharpen  the  Doppler  estimates.  Also,  the  absolute 
Doppler  information  is  needed  in  the  fast  track  solver  described  in  Appendix 
3. 
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A. 11  Fast  Estimation  of  Track  Parameters 
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A  further  computational  advantage  can  be  obtained  by  using  a  path 
evaluation  method  which  is  less  comprehensive  than  BEARTRK  when  the  goal  is  to  y- 

solve  the  path  association  problem  across  time  and  across  sensor  pairs.  Once  S’* 

the  association  problem  is  solved,  BEARTRK  can  be  called  to  compute  optimum 

JP 

track  estimates  from  single-source  path  data  for  each  source.  A  fast  track 
solver  is  presented  in  Appendix  B. 
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Appendix  B  -  Fast  Estimation  of  Track  Parameters 

When  tracking  multiple  sources,  it  may  be  necessary  to  try  many  different 
associations  of  spectral  peaks  to  lines  and  sources.  In  such  a  case,  the 
ability  to  quickly  estimate  accurate  track  parameters  is  essential. 

B.l  A  Fast  Track  Solver  for  the  Single-source,  Single-Velocity  Case 

The  algorithm  is  a  modified  version  of  the  one  described  in  the  SCT  85- 
371  proposal.  The  principal  modification,  suggested  by  R.  Bliss,  is  the 
removal  of  the  assumption  that  the  range  at  CPA  be  comparable  among  the 
various  sensors. 

The  algorithm  works  by  measuring  "S-curve"  parameters  and  inferring  track 
parameters.  Ry  “S-curve"  we  mean  the  general  appearance  of  the  Doppler  vs. 
time  observed  at  a  single  sensor  when  a  source  is  passing  by  at  a  fixed  speed. 

In  order  to  use  the  algorithm  with  the  OTS  system,  it  is  necessary  to 
convert  from  intersensor  Doppler  (as  measured  by  CORAN,  for  example)  to  an 
estimate  of  absolute  Doppler.  Alternatively,  an  "auto-Doppler  correlogram  " 
(ADC)  can  be  computed  directly,  wherein  all  spectrograms  for  a  given  sensor 
are  Doppler-correlated  against  the  spectrogram  at  CPA  in  that  sensor;  this 
will  produce  what  we  call  an  S-curve.  The  CPA  frequencies  measurable  in  the 
spectrogram  can  be  used  to  convert  the  "cross-Doppler  correlograms"  (CDC) 
produced  by  CORAN  into  equivalent  ADC's,  thus  providing  ADC's  with  the  higher 
noise  immunity  inherent  in  cross-correlations. 

The  output  of  ADEC  (frequency  tracks  versus  time)  can  be  used  directly  to 
provide  S-curves.  Dividing  each  line  by  its  CPA  frequency  and  averaging  those 
which  appear  to  coincide  (because  they  are  from  the  same  source)  can  be  used 
to  increase  noise  immunity.  Yet  another  alternative  is  to  produce  absolute 
Doppler  estimates  using  MAPLE. 

8.2  Algorithm  Description 


Consider  the  case  of  Nr  receivers  and  N$=l  sources.  The  source  is 
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J 


assumed  to  be  heading  on  a  fixed  course  and  speed  without  acceleration.  Given  '• 

the  ADC's,  there  will  exist  exactly  one  S-curve  for  each  sensor  as  shown  in 

<-i( 

Fig.  1.  (In  the  multisource  case,  an  isolated  S-curve  is  hypothesized  by  some  .y 
choice  of  correlogram  peak  association  over  time.  In  the  case  of  a  single 

r." 

source  with  multiple  sinusoidal  lines,  the  S-curves  can  be  combined  by  r.| 

averaging. ) 
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Figure  1.  Doppler  Shift  vs.  Time  for  a  Constant-Velocity  Source 

The  S-curve  is  expressed  as 

d(t)  =  1  -  ^  cos[9y(t)  -9-(t)]  =  Doppler-shi ft  vs.  time  (1) 

where  v  is  source  speed,  c  is  sound  speed,  a .  is  the  bearing  of  the  source 
from  sensor  i,  and  0  is  the  bearing  of  the  source  track,  as  shown  in  Fig. 

2. 
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It  is  assumed  that  the  sensor  positions  and  the  speed  of  sound  are  known 
constants.  From  the  S-curve  in  sensor  i  (Fig.  1)  we  compute 


d+  =  max  d. (t) 
t 

dT  =  min  d.  (t) 
1  .  1 


t^  =  Time  of  CPA  at  sensor  i  =>  d ( t ^  )  =  1 


d+'$ 


3d  (t) 

i  1 


at 


■( t  ^ )  |  =  max 


3d  (t) 

max|~^ - (s) 


(2) 


Inspection  of  (1)  reveals  that  for  a  source  speed  of  v  the  Doppler  shift  is 
between  d  =  (1  +  v/c)  and  d  =  (1  -  v/c)  .  Thus,  the  speed  of  the  source 
is  estimated  from  each  S-curve  as 


d  -d. 

v  =  c  J _ 1 

1 

di  *  di 


and  the  source  speed  estimate  is  given  by  a  simple  average: 


(3) 


Equation  (1)  also  shows  that,  at  CPA, 


di  =  cTT 


R 


Ri 


-  1  v 

"  TT  l 

R  1=1  cd. 


From  Fig.  3,  we  find  that  the  angle  between  the  track  and  the 
between  sensors  i  and  j  can  be  computed  from 


s,n<Yu>  ■ 


1J 


(5) 


where  is  the  source-recei ver  range  at  CPA  for  the  i-th  receiver, 
estimate  of  range  at  CPA  in  each  sensor  is  obtained  as 

N„ 


Thus,  an 


(6) 


i  ne 


(7) 


where  R^j  is  the  distance  between  sensors  i  and  j.  Note  that  from  these 
measurements  the  sign  of  the  angle  y„  is  ambiguous.  If  denotes  the 

angle  of  the  line  from  sensor  i  to  sensor  j,  then  the  angle  of  the  track  can 
be  estimated  as 


ev(i,j)  *  t  ^ j 


Figure  3.  Two-Sensor  Geometry 


Averaging  over  all  pairings,  and  choosing  a  consistent  set  of  signs  in  (8),  we 
obtain  the  track  angle  estimate  as 

A  A  1  A 

ev  =  (N„  -  1J!  I  ^ 

ij 

A 

The  values  of  ey ( i , j )  can  be  examined  for  outliners  and  pruned  accordingly. 


From  the  measured  CPA  times  ti ,  and  the  estimated  ranges  at  CPA  Ri ,  we 

/*. 

can  use  the  track  angle  estimate  to  find  an  estimate  R  (i)  of  the  range  R0 

of  the  source  from  the  origin  (0,0)  at  the  time  of  CPA  t  with  respect  to 

0 

the  origin.  The  (polar)  coordinate  system  for  this  is  shown  in  Fig.  4;  in 
these  coordinates,  the  track  is  the  line  orthogonal  to  the  vector  from  the 
origin  to  the  point  (R  ev  +  it/2)  .  Let  sensor  i  have  absolute  polar 
coordinates  (r  ,  ^ )  .  Then 

PqO')  «  -  riSin(Vh)  (10) 


The  average  over  all  sensors 
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gives  an  estimate  of  the  range  from  zero  (0,0)  at  CPA. 
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Figure  4. 


Absolute  °olar  Coordinates 


Finally,  given  the  angle  of  the  track  and  the  range  Rq  at  CPA  with 
respect  to  the  origin,  the  time  of  CPA  to  the  origin  tq  can  be  estimated  by 
averaging  the  CPA  times  t ^  at  each  sensor  minus  the  time  of  travel  from  CPA  at 
the  origin  to  CPA  at  sensor  i: 

1  Nr  ri  C0SOi  -  0„) 

**  - -  <'?> 

R  1=1 


The  resulting  parameters  (Rq,  e  T  ]  are  the  track  parameters  output 
by  the  algorithm. 


Appendix  C  -  Measurements  for  Localization 


In  this  appendix  we  describe  measurements  which  are  useful  for  multi  - 
source  localization.  The  original  data  consist  of  NR  digitized  acoustic 
hydrophone  signals. 

C.l  Primary  Measurements;  Auto-  and  Cross-Spectra 


All  information  is  derived  from  the  measured  power  spectral  density  (PSD) 
S-  t{u)  in  each  sensor  and  the  cross-spectral  density  (CSD)  S.  .  (u) 

between  each  pair  of  sensors.  The  PSD  is  computed  as  described  in  [1],  and 
the  CSD  computation  is  exactly  analogous  [2].  Briefly  summarized,  the  PSD 
S-  (u>)  is  an  average  of  squared-magnitude  FFT's  of  successive  blocks  of 
digitized  pressure  versus  time  in  sensor  i.  Similarly,  the  CSD  is  an  average 
of  the  FFT  at  sensor  i  times  the  conjugate  of  the  FFT  at  sensor  j.  Frequency 
u  =  2irf  is  in  radians  per  second  and  is  regarded  as  continuous  via 
interpolation  of  FFT  bins  when  necessary.  Time  t  is  in  seconds  and  is  defined 
as  the  mid-point  of  the  averaging  interval  used  to  measure  the  PSD  or  CSD. 

Time  is  also  available  continuously  through  interpolation. 


The  time-varying  PSD  S.  (u)  ,  regarded  as  a  plot  of  power  density 

1  yZ 

versus  time  and  frequency,  with  i  fixed,  will  be  called  the  spectrogram  of  the 
data  from  sensor  i.  Similarly,  S...  fc(u>)  will  be  called  the  (complex)  cross- 

*  j ,  t 

spectrogram. 


C.2  Secondary  Measurements:  Auto/Cross  Delay/Doppler  Correiograms 

Recall  that  the  auto(cross)-correlation  function  is  equal  to  the  inverse 
Fourier  transform  of  the  power(cross)-spectral  density.  Accordingly,  we 
define  the  delay  auto(cross)  correl ogram  as  the  two-dimensional  distribution 
obtained  by  replacing  each  slice  along  frequency  (time  fixed)  of  the  cross- 
spectrogram  by  its  inverse  FFT.  The  term  "correlogram  "  alone  will  refer  to  a 
delay  cross-correl ogram.  The  auto-correl ogram  can  be  useful  for  tracking 
multipath  delays. 


C.3  Differentia  1  Delay  Correlograms 


The  correlogram  between  sensors  i  and  j  based  on  data  centered  at  time  t 
is  given  by 

-  i  if  jut  . 

-TT 

where  t  is  the  correlation  lag.  The  maximum  of  r--  .  f 1 1  with  respect  to 
-  i  J  >  t 

t  is  normally  defined  as  the  time-difference  of  arrival  (TDOA)  between 

sensors  i  and  j. 

C.4  Differential  Doppler  Correlograms 


The  second  derived  distribution  is  called  the  Floppier  correlogram.  The 
Doppler  auto(cross)-correlogram  between  sensors  i  and  j  based  on  data  centered 
on  time  t  is  defined  by 

DU,t(6)  4  t  (13> 

“TT 

where  3  is  called  the  differential  Doppler  between  sensors  i  and  j  at  time 

t.  We  have  D..  .(3)  =  D..  .  ( 1 /0 ) /s  .  The  Doppler  auto-correl ogram  can  be 
1  j  ,  1-  J  i ,  t 

denoted  more  simply  by  0.  ^(3)  =  D. .  ^(3)  . 

C.5  Computation  of  the  Doppler  Correlogram 

The  Doppler  correlogram  can  he  efficiently  computed  using  an  FFT 
correlation  facility  [3].  The  basic  idea  is  to  sample  the  power  spectra 
versus  log  frequency  in  order  to  convert  the  "stretch"  operation  applied  to 
S-  t(u>)  in  (13)  to  a  simple  shift  operation.  This  converts  equation  (13) 
into  a  normal  cross-correlation. 


The  power  spectra  in  equation  (13)  are  transformed  by  the  change  of 
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to  yield 

\«H4si,.W-si/l 

Equation  (13)  becomes 


0u>t(8)  ■  /  *»>*“« 


which  is  simply  the  cross-correlation  of  the  exponentially  sampled  power 
spectra  in  channels  i  and  j,  weighted  by  frequency  w  =  exp(w)  .  The  weight 
factor  can  be  absorbed  into  the  power  spectra  so  that  an  unweighted  cross- 
correlation  facility  can  be  used: 

n..  (b)  -JL  j11  s'  tU)s;  (b  +«)«:  (17) 

/?  -TT  * 


where 


Sj.t®  s  «“/zJ,.tOT  -- 


C.6  Spectral  Resampling 

In  practice,  the  power  spectra  are  available  in  sampled  form 
Si  t^k)  =  2*k/N’  k=fM»2,...,Nf-l  ,  where  Nf  is  the  number  of  uniformly 
spaced  samples  around  the  unit  circle  (typically  the  FFT  size).  In  this 
situation,  the  exponentially  sampled  replacement  T.  t(uk)  is  computed  using 
digital  interpolation. 

Typically,  the  power  spectrum  is  assumed  to  be  bandl imited  to  less  than 
one-half  the  sampling  frequency  fs.  Unfortunately,  this  requires  that  the 
corresponding  autocorrelation  function  be  of  infinite  extent  in  the  lag 
domain.  Spectral  interpolation  is  meaningful  only  when  the  corresponding 
correlation  function  asymptotically  approaches  zero  (excluding  periodic 
correlation  functions  which  are  easily  interpolated  with  zeros  between  the 
original  spectral  samples).  Therefore,  proper  spectral  interpolation  is 
carried  out  by  first  extending  the  autocorrelation  function  by  means  of 
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bandlimited  extrapolationm  [4]  until  the  correlation  function  falls  below  the 
minimum  representable  number  (i.e.,  it  is  quantized  to  zero  beyond  some  lag), 
and  then  performing  time-limited  interpolation  [5]  in  the  frequency  domain. 
Time-limited  interpolation  in  the  frequency  domain  is  accomplished  by  zero¬ 
padding  the  extrapolated  correlation  function  and  taking  the  FFT. 

A  method  for  power  spectrum  interpolation  which  avoids  bandlimited 
extrapolation  is  as  follows.  Assume  the  autocorrelation  function  is  Bartlett 

windowed.  That  is,  if  r.  fc(r)  denotes  the  sample  autocorrelation  function 

——————  i , t 

for  sensor  i  at  time  t,  the  Bartlett  (or  triangular)  windowed  autocorrelation 
is  given  by 


Nf-2|- 


t(xi  «ir,-t(T)  Nf 

1  » u  n 


,  |t|  <  Nf/2 

,  | t I  >  Nf/2 


The  power  spectrum  then  is  the  convolution  of  the  true  power  spectral  density 
(PSD)  with  a  sine2  function  (defined  below)  having  a  zero-crossing  interval 
equal  to  4  spectral  samples  (or  "bins"). 

When  the  periodogram  method  for  PSD  estimation  is  used  [7],  as  we  use 

here  [1],  a  length  Nf/2  rectangular  window  is  applied  to  the  data  segment,  a 

length  Nf  FFT  is  taken  (with  zero  padding),  and  the  magnitude  squared  is 

averaged  over  successive  blocks.  In  this  case,  the  true  PSD  is  convolved  with 
2 

the  same  sine  function  as  when  Bartlett-windowing  the  true  correlation 
function . 

Because  the  periodogram  method  of  power  spectrum  estimation  is  an  average 
of  the  squared  magnitude  of  finite-length  transforms,  time-limited 
interpolation  is  a  valid  operation  (assuming  sufficient  zero-padding  so  that 
the  operation  of  taking  magnitude  squared  of  the  FFT  results  in  a  non-circular 
Bartlett-windowed  sample  autocorrelation). 

Given  periodogram-method  PSD  samples  S..  (u^),  =  2*1 /N^.  , 

k  =  0,1,2,...,  N^-l,  the  interpolated  values  are  defined  by 


\  \  \  v 
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si,t {a)  si,t(wk)sinc(“  ‘  “k> 


where 


sinc(o))  - 


sin{uj/2)Nf 


In  practice,  the  sine  interpolating  function  is  normally  symmetrically 


windowed  to  the  first  K  -  3  zero  crossings  in  each  direction.  If  w  (»,  ) 

s  s  k 

denotes  this  window  function  of  length  2K$  +  1  spectral  bins,  then  the 


interpolated  values  are  approximated  by 


Si  t^  “  l  Si  t )<*>$ (“  "  “|c)sincl(»  "  \) 
*  k=k^  ’ 


where 


uj  (uj)  =  0.54  +  0.46  cos  (w  )  ,  u  <• 


(Hamming  window) 


kj  -  max {0, 


<2  =  min{Nf/2,  u»  +  u>$: 


“s  =  2lTVNf 


«s  =  Number  of  zero-crossing  retained  per  wing 


If  the  intersample  spacing  is  to  be  greater  than  the  original  sample 
spacing,  it  is  desirable  to  smooth  the  power  spectrum  so  that  aliasing  does 
not  occur  in  the  lag  domain.  Let  T '  >  1  denote  the  desired  maximum  bin 
spacing  in  the  frequency  domain.  (  T '  =  I  yields  the  original  bin 
spacing.)  The  ideal  interpolation  formula  (20)  extended  to  allow  subsampling 
without  aliasing  is  then 


Si,t^  "  T'  Si  ,t^k^sinC'-  T1 


The  practical  windowed  version  is  analogous  to  (22).  An  efficient  table- 
lookup  implementation  of  windowed-sine  interpolation  is  described  in  [6]. 


C.7  Exponentially  Spaced  Spectral  Sampling 


We  prepare  the  sample  PSD  versus  sampled  log  frequency  TT.  (u)  by 

T  >  t 

resampling  the  linearly  spaced  frequency  axis  to  produce 

It  is  necessary  to  choose  the  total  frequency  interval  and  sampling  density 
over  which  to  perform  this  transformation.'  For  convenience  in  later 
processing,  we  specify  (1)  the  highest  frequency  to  be  retained  after 
mapping  and  (2)  the  number  TT^.  of  resulting  samples  in  TT  (^)  .  The  lower 
cutoff  frequency  is  computed  from  these  inputs  along  with  the  maximum  allowed 
bin  spacing  which  is  always  set  to  the  original  bin  spacing  (to  avoid  loss  of 
spectral  information  in  the  retained  band). 

The  resulting  algorithm  is  as  follows.  Let  P(k)  denote  the  array 
containing  Si  tUk),  k»0,l,...,Kf  <  Nf/2  ,  and  let 


y(k)  *  P(Frk), 


TT=0, 1 , . . .  ,77- 1 


denote  the  corresponding  resampled  array,  where  kQ  is  the  bin  coordinate  of 
the  lowest  frequency  retained,  and  r>l  is  the  ratio  of  the  sampling  interval 
increase.  The  requirement  ^  ^  <  2tt /N^-  reduces  approximately  to  the 

requirement  ui^(l-l/r)  <  2*/Nj.  which  in  conjunction  with  specification  of 
and  TT  yields 

k, 

.  A  ,  -TT 

ko  *  V 


where  k^  <  N^/2  is  the  bin  coordinate  of  the  highest  frequency  retained. 
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OF  SOURCE  LOCALIZATION 

J.  0.  Smith  and  J.  S.  Abel 

Systems  Control  Technology 
1801  Page  Mill  Road 
Palo  Alto,  CA  94304 

ABSTRACT 

A  closed-form  least-squares  approximate  maximum  likelihood  method  for 
localization  of  broad-band  emitters  from  time-di f ference-of-arri val  (TDOA) 
measurements,  called  the  Spherical  Interpoal ti on  (SI)  method,  is  presented, 
i he  localization  formula  is  derived  from  least-squares  "equation-error" 
minimization.  Computer  simulation  results  show  that  the  SI  method  has 
variance  approaching  the  Cramer-Rao  lower  DOund. 

1.  Introduction 

The  problem  of  automatically  locating  a  radiating  or  reflecting  source  by 
analyzing  signals  received  from  the  object  is  a  basic  one  having  numerous 
applications  in  defense,  aerospace,  geophysics,  and  industry.  This  paper 
discusses  a  highly  efficient  solution  to  the  problem  of  locating  a  source  with 
a  passive,  stationary  sensor  array. 

A  source  reveals  informat. an  about  its  location  througn  the  relative  time 
delays  seen  among  the  signals  it  radiates  to  a  sensor  array.  In  a  constant- 
velocity  medium,  the  time  difference  of  arrival  (TnOA)  between  signals 
received  ’ n  a  sensor  pair  will  place  the  source  on  a  hyperboloid  of  revolution 
with  an  axis  along  tne  line  drawn  between  *.ne  sensors.  In  an  n-cimensiona1. 
space,  n  suer,  "nOA;s  fron  n  nondegenerate  1  y  placed  sensor  pairs  are 
necessany  and  sufficient  to  uniquely  determine  tne  source  location  [si. 

Source  localization  can  therefore  be  performed  using  n-i  sensors  placed  so 
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that  they  do  not  lie  in  a  subspace  of  dimension  less  than  n.  However,  to 
miminize  the  effects  of  noise  in  the  TDOA  measurements,  it  is  desirable  to  use 
N  >  n+1  sensors  and  fit  to  all  of  the  TDOA's  by  minimizing  some  measure  of 
inconsistency  between  the  measured  TDOA's  and  the  TDOA's  implied  by  the 
estimated  source  location. 

While  an  extensive  literature  exists  on  the  problem  of  estimating  TDOA's 
from  received  signals  [4],  very  few  papers  seem  to  be  available  on  the  problem 
of  converting  TDOA  values  into  source  location  [1,2,3,5,11,12]. 

Conventionally,  source  locations  are  estimated  by  intersecting  hyperbolic 
lines  of  position  (LOP)  determined  by  range  difference  measurements 
[5,11,12].  However,  finding  the  intersection  of  a  set  of  hyperboloids  is 
computationally  intensive,  involving  finding  the  minimum  of  a  non-convex 
function.  R.  0.  Schmidt  [1]  has  proposed  a  formulation  in  which  the  source 
location  is  found  as  a  focus  of  a  conic  passing  through  three  sensors. 
Schmidt's  method  can  be  extended  to  an  optimal  closed-form  localization 
techniaue  [8].  In  [2],  a  gradient  search  localization  procedure  is  derived 
for  computing  optimal  source  locations  from  noisy  TDOA's.  In  [3],  a  formula 
is  given  for  single-source  TDOA  localization  from  four  sensors,  and  can  oe 
extended  to  an  arbitrary  number  of  sensors  [8]. 

In  [8],  a  closed-form  localization  technique,  termed  the  Spherical  - 
Interpol ati on  (Si)  method,  was  described  and  was  shown  to  perform  better  than 
two  related  techniques  [1,3],  In  this  paper,  we  present  the  SI  method  and 
develop  expressions  for  the  variance  of  SI  source  location  estimates  in  the 
presence  of  noisy  TDOA  measurements.  We  give  a  geometric  interpretation  of 
the  source  location  estimates  produced  by  the  SI  method  and  show  that  the  SI 
estimates  are  related  to  maximum  likelihood  estimates.  In  addition, 
simulation  results  are  presented  in  which  the  SI  method  shows  noise  immunity 
approaching  the  Cramer-Rao  lower  bound. 

The  structure  of  the  paper  is  as  follows.  In  section  2  we  derive  tne  SI 
method  for  closed-form  localization  of  a  source  in  a  field  of  ‘1  sensors.  In 
section  3,  simulation  results  are  presented  for  two  different  source  locations 
and  two  additive  TDOA  noise  levels.  Section  a  gives  a  geometric 
; nterpretation  of  the  SI  technique. 
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2.  Closed-Form  Localization  of  One  Source  Using  N  Sensors 


Let  N  denote  the  number  of  sensors,  and  let  d. ,/c  denote  the  TDOA 

ij 

between  sensors  i  and  j  (  i,j=l,...,N  )  where  c  is  the  speed  of 

propagation.  The  vector  of  (x,y,z)  spatial  coordinates  for  the  ith  sensor 

is  denoted  x.,  and  the  position  of  the  source  is  denoted  x  .  The  distance 
— i  — s  - 

between  the  source  and  sensor  i  is  denoted  by  D.  =  »x.  -  x  a,  and  the 

i  i  — s 

distance  from  the  origin  to  the  point  x.  is  denoted  R  .  Similarly. 

—i  i 

Rs  *  jx^il.  These  quantities  appear  in  Figure  la. 

The  TDOA  between  sensors  i  and  j  is  equal  to  the  range  difference 
(RD)  d^  divided  by  the  speed  of  propagation,  a  constant  velocity  medium  is 
assumed.  It  will  be  more  convenient  to  work  directly  with  RD's  instead  of 
TDOA' s. 


The  localization  problem  is  to  determine  x  given  d  for  i  and  j 

ij 

between  1  and  N.  flote  that  there  are 

rn,  N(N-l) 

2 J  2 

distinct  Rfl's  d..  (excluding  i=j,  and  counting  each  d..  =  -d . .  pair  once) 
ij  ij  Ji 

however,  any  N-l  RD  measurements  which  form  a  "minimal  spanning  subtree" 

determine  all  the  rest  (in  the  noiseless  case).  The  redundancy  of  the 

complete  set  of  RD  measurements  is  used  to  increase  noise  immunity. 

We  have  the  following  basic  relations: 


dij  =°i  -°r 


i  =  1, . . .  ,n,  j=i ,. . . ,n 


°i  *  !'^i  *  —5 


0?  =  I X  .  1 2  -  2<x .  ,x_>  +  IX 


2.1  The  Equation-Error  Formulation 


We  assume  here  that  the  set  of  RD  measurements  d  ,  i=l,...,N,  for  some 
j,  are  available.  To  solve  the  localization  problem,  we  first  map  the  spatial 
origin  to  an  arbitrary  sensor,  say  the  jth.  This  gives 


x.  =  4  n 

-J 


I V  " 

)  °j  •  8s 


From  (1)  and  (4)  we  have 


0.  =n.*d..=R  +  d . .  (5) 

i  J  IJ  s  ij 

as  illustrated  in  Figure  lb.  Substituting  R^  +  d.  .  for  n.  in  ( .1)  yields 

R2  -  d2.  -  2R  d.  .  -  2xTx_  =  C:  (n) 

i  i j  s  ij 

The  jth  equation  is  degenerate  so  we  have  N-l  eauations  in  three  unknowns 


As  the  delays  are  typically  not  measured  precisely,  we  introduce  a  so- 
called  "equation  error"  [7]  into  the  ri ght-hand-si de  of  (fi),  and  minimize  it 
in  a  least  squares  sense  to  provide  an  estimate  of  the  true  solution.  Without 
loss  of  generality,  let  j=l.  Then  (A)  becomes 


=  R2  -  d2  -  2R  d . ,  -  2xTx  , 
i  i  U  s  il  — i— s  ’ 


i  =  2,3,...  ,M 


wnere  e.  is  the  equation  error  to  be  minimized  The  set  of  '!-!  equations 
(7)  can  be  written  in  matrix  notation  as 


e»3-2Rd-2Sx 
-  -  s-  -s 


i 


It  is  worth  noting  that,  formally,  equation  (8)  is  linear  in  x  given  R  , 


and  it  is  linear  in  R^  given  Error  vectors  which  are  linear  in  the 


unknowns  yield  closed-form  least-squares  solutions. 


The  formal  least-squares  solution  for  x  given  R  is 

— s  s 


is  '  J  S*(i  -  2S,aj 


(?) 


where 


s*  =  f  $Ts  I'V 


(in) 


yields  the  unweighted  least-squares  solution.  If  it  is  desired  to  weight  the 


RD's  according  to  a  priori  confidence  in  each  RD,  then  the  weighted  equation 


error  energy  e  W  z  is  minimized  for 


*  T  -IT 

s  -  c  s'w  S  )  s  * 


(11) 


where  W  is  positive  definite. 


To  obtain  a  true  least  squares  estimate  for  ,  it  is  necessary  to 

minimize  j  =  (or  j  =  G'w  e)  with  respect  to  while  allowing  R^ 

to  vary,  maintaining  the  relation  R  =  ;ix  i . 

$  — s 

Unfortunately ,  this  minimization  is  over  a  ron-c onvex  cost  'unction,  anc 
exnaustive  search  techniques  must  be  used  in  general.  However,  in  this  case 
closed-form  solution  can  be  found  wnioh  aoDroximately  minimizes  J  =  . 


2.2  The  Spherical-Interpolation  Method 


The  basic  idea  of  the  new  closed-form  solution  is  to  substitute  (9)  into 

(8)  and  minimize  the  equation  error  again,  this  time  with  respect  to  R  . 

This,  surprisingly,  yields  a  linear  least  squares  problem  for  finding  R  ,  an 

the  solution  is  computationally  inexpensive.  The  technique  is  made  possible 

by  the  fact  that  the  formal  least  squares  estimate  of  given  Rs  in  (8) 

is  itself  linear  in  R  .  When  the  minimizing  R  value  is  found  in  this  new 

s  s 

linear  equation,  the  corresponding  value  of  x  (via  (9))  is  automatically  a 

— s 

minimizer  of  the  squared  equation-error  norm  with  respect  to  x  given  this 

— 5 

V 


Rewriting  the  equation  error  (8)  to  eliminate  x  by  subsituting  the 

— s 

value  from  (9),  we  get 

e  -  i  -  29^  -  S  S*  (5.  -  2Rsd]  «  ( I  -  S  S*  ) (a  - 

where  I  is  the  N-l  by  N-l  identity  matrix.  Now,  £  is  1  inear  in  the 
single  unknown  R  .  nefine  the  N-l  by  N-l  symmetric  matrix 

a  *  T  -IT 

Ps  =  S  S  =  S  (  S  W  S  )  S  M 

The  rank  of  is  at  most  3  regardless  of  its  order  N-l.  Also,  P,.  is 
2 

idempotent  (  p$  s  Ps  )  *  p i ^ 3 1 1 y ,  the  identity  matrix  minus  an  idempotent 

9 

matrix  is  idempotent,  i.e.,  (  I  -  =  I  -  P^  .  Idempotent  matrices  are 

projection  operators;  the  operation  P$  x_  will  remove  components  of  _x  n°T 
lying  in  the  space  spanned  by  the  columns  of  S. 

In  the  four-sensor  case,  P^  =  I,  and  the  error  £  is  the  zero  vector. 

In  the  more  general  case  of  N  >  4  sensors, 

e  *  Ps  U  -  2Rsd)  =  (I  -  P-  )(i  -2Rsdj  (12) 

so  that  the  equat i on-error  energy  becomes 


J  3  £L  s  (1  ‘  2Rci)T  T  («.  -  2R  d) 


where, 


T  =  PX 
1  KS 

(or,  to  minimize  ,j  =  eTW  e_  ,  let  T  ^  P^W  P^  )  .  Minimizing  J  with  respect  to 
R  is  a  form  of  weighted  least  squares  in  which  the  weighting  matrix 
T  is  idempotent  of  rank  N-4.  The  missing  dimensions  reflect  the  degrees  of 
freedom  removed  by  choosing  sensor  1  as  the  origin  and  substituting  in  the 
least-squares  solution  (9)  for  the  three  spatial  source  coordinates.  The 
solution  is  given  by 

,  d T  T  a 

Rs  =  7  — —  d3) 

P.d T  d 

Substituting  this  solution  into  (9)  yields  source  location  estimate 

is  *  \  s*u  -  2Rsi)  d4) 

Clearly,  the  computati onal  burden  of  (14)  is  very  low  compared  to  iterative 
nonlinear  minimization.  If  iterative  nonlinear  minimization  is  desired  (to 
obtain  the  lowest  possible  variance  and  bias),  (14)  provides  an  excellent 
initial  value  for  a  general  descent  method. 

a  • 

note  that  the  pair  (x  ,  R  )  minimizes  the  equation-error  energy 
t  — s  s  .  . 

j  =  e'  w  e  without  the  constraint  s  ax  :i  •  We  expect  that  the  pair 


,  R$  =  ii x ^ ii  1  approximately  minimizes  J  subject  to  the  constraint 
jx  j  *  R  .  Therefore,  we  define  the  range  estimate  by 


Rs  =  ix,, 


instead  of  using  (13).  Similarly,  the  bearing  estimate  is  defined  as  the 
vector  of  direction  cosines  from  the  origin  at  sensor  1  to  the  source: 


x  ] 

-s 


V2 


Note  that  the  SI  solution  (14)  is  based  on  RDs  measured  relative  to  a 
single  reference  sensor.  When  additional  RDs  are  available,  additional  SI 
source  location  estimates  can  be  made,  using  different  reference  sensors,  and 
the  source  location  estimate  could  be  computed  as  the  weighted  average  of  the 
SI  estimates. 


2.3  Variance  of  the  SI  Estimate 


Here,  we  give  expressions  for  the  variance  of  the  unweighted  SI  estimates 
of  source  location,  range  and  bearing.  Assuming  the  RD  variance  is  small 
compared  to  the  RD  mean,  the  variance  of  the  SI  estimator  of  source  location 
is  given  by  [10] 


''■“-(is)  ■  Cj?)  «„(#) 


where  the  RD  vector,  d  is  assumed  unbiased  with  covariance  matrix  R 


Accordingly,  the  variance  of  the  source  range  and  source  direction  cosines 
vector  are  given  by 


A  3R  „  3R 

varCRs)  '  lj£)  '’•'•(i.HyV 


a  a  t 

var(&)  -  !#)  wrfx  )(t?)t 


'•-s;''3x  ■ 

— s 


The  derivative  3*  /3d.  can  be  evaluated  by  using  (8)  in  the  noiseless  case, 


=  0: 


3_6  3R  _d  3x 

3d  •  2  ST  -  2STd  =  0 


Recall  ^  =  R*  -  d^  >  and 


17  =  _2Ad 


wnere 


«. «. ' 


-•  u"  v 


The  derivative  may  be  evaluated  using  the  chain  rule: 


3R-d  aR, 

5—  _  Af  5 


T  /  S 


~W  =  +  Rs  a?  =  -  Ss^)  +  V 


where  R$  £  iix_s<j  and  £$/Rs  *  £s  have  been  used  to  evaluate  '3Rs/3x_s  . 
(19)  and  (20)  in  (18)  we  find. 


Using 


where 


o  X 

Ad  *  Rs 1  +  Cl  Us  +  S  bar  = 


ir-  -(  )"1aT(  *d  +  V  ) 


A  -  (d.  S  ) 


Using  the  definitions  for  ^  and  ,  we  have, 


air;  1  — s 


-s  1  _i  il,.  T, 

3s£  =  R7  Pxs  =  1  -  ilsiis ) 

where  P*  represents  a  orojection  operator  which  removes  components  of  a 
-s 

vector  along  the  direction  of  x  ,  Substituting  equations  (21)  and  the  above 

— $ 

into  (16)  and  (17),  we  have  the  desired  result: 


Var(iLs)  '  UTa  )aT  (*j  +  R$I  )^d( Ad-RsI  jA  f  aTa 


7ar(Rsj  -  ^  Varrx.)^ 


VaPCSsj  '  h  V  Var(*s)px 

s  - s  ~s 


3.  Simulation  Results 


This  section  presents  simulation  results  on  the  performance  of  the  SI 
method  on  the  problem  of  localizing  a  single  source  in  the  presence  of 
noise.  The  simulations  were  implemented  in  the  Ctrl -C*  language  on  a  Vax 
11/780  computer.  The  performance  of  the  SI  localizer  is  evaluated  according 
to  the  bias  and  variance  of  the  source  location  estimate.  The  results  show 
that  the  variance  of  the  SI  method  is  given  accurately  by  (22),  and  that 
performance  of  the  SI  method  is  nearly  optimal,  i.e.,  the  expected  square 
error  in  the  source  location  estimate  approaches  the  Cramer-Rao  lower  bound. 

All  simulations  employed  the  9-sensor  array  geometry  shown  in  Figure  2, 
with  the  two  source  locations  shown  and  two  levels  of  additive  white  noise  in 
the  RD's.  The  sample  bias  and  variance  were  obtained  by  averaging  the  results 
of  1000-trial  Monte-Carlo  runs. 

Typically,  RD  estimates  are  zero  mean  with  variance  dependent  only  on  the 
signal-to-noise  ratio  (SNR)  at  each  sensor  [A, 6].  When  the  source-to-sensor 
range  is  comparable  in  all  sensors,  the  SNR  is  normally  comparable  also.  Our 
simulated  range  difference  estimates  were  generated  by  adding  white  Gaussian 
noise  to  true  RO  values,  corresponding  to  the  case  of  uniform  SNR  across 
sensors. 


I  ' 

*  *  Ctrl-C  is  a  registered  trademark  of  SCT. 
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Table  1.  Envi romental  information  for  the  four  simulation  cases. 
(Note:  all  distances  are  measured  in  meters.) 


Run  1 


Run  2 


Source  location  (390,160,170) 

Range  Rs:  454.5 

Bearing  cosines  a  :  (.858, .352, .374) 

—-5 

R0  noise  standard  deviation  g.  =0.1 

du 

Source  location  (390,160,170) 

Range  R$:  454.5 

Bearing  cosines  n  :  ( .858, . 352, . 374) 

“S 

RO  noise  standard  deviation  g^  =  1.0 


Run  3 


Run  4 


Source  location  _x  :  (540,1360,110) 

Range  R^ :  14/57 

Bearing  cosines  2  '■  (•  368,  .827 ,  .075) 

RD  noise  standard  deviation  jd  =  0.1 

Source  location  x  :  (540,1360,110) 

Range  Rs:  1467 

Bearing  cosines  i  :  (.368, .927, .075) 

—"S 

RO  noise  standard  deviation  g.  =  1.0 

du 


Table  1  describes  the  enviromental  information  (source  location,  range, 
bearing,  and  additive  noise  level)  for  each  of  four  1000-trial  Monte-Carlo 
runs  (two  source  locations  and  two  TOOA  noise  levels).  For  all  runs,  the  nine 
sensors  were  located  at  (0,0,0),  (0,0,100),  (0,0,200),  (100,0,0),  (100,0,100), 
(100,0,200),  (0,100,0),  (0,100,100),  and  (0,100,200)  meters. 
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Table  2.  Sample  bias  measurements,  SI  Method  1000  trials 


■  F 


Table  3.  Sample  standard  deviation  measurements,  SI  Method  1000  trials 


!>■ 

f' 

Run 

Source  Location 

Range 

Source  Bearing 

(xl0d) 

at 

°y 

Js 

°z 

s 

°R 

s 

y 

Sample  1 

1.9 

0.66 

0.42 

2.0 

5.6 

7.0 

9.5 

Standard  2 

19. 

6.5 

4.1 

20. 

56 

70 

95 

Deviation  3 

9.5 

26. 

0.58 

27. 

7.3 

6.4 

13 

' 

\- 

4 

32. 

220. 

5.2 

235. 

78 

63. 

126 

Standard  1 

1.7 

0.62 

0.30 

1.3 

4.5 

5.0 

8.7 

£ 

Deviation  2 

17. 

6.2 

3.0 

18. 

45. 

50. 

87. 

from  (22)  3 

9.3 

25. 

0.44 

26. 

6.2 

3.0 

12. 

r 

4 

93. 

248. 

4.4 

264. 

62. 

30. 

123. 

k 

Table  3  snows  the 

sample  standard  deviations  for  eacn  Monte-Carlo  run, 

again  for  both  Euclidean  and  plane 

-projected  polar 

coordinates  and 

range.  Due 

to  the  source-sensor  geometry,  the 

bearing 

is  more 

accurately 

estimated  than 

the  range  in  all 

examples;  this  is 

typical 

when  the 

source  is 

several  aperture 

sizes  away  from 

the  sensor  array. 

Note  the 

estimate  variance 

of  th 

e  SI  method 

,V 

appears  to  increase  linearly  with  , 

3n  increase  in  RD 

vari ance. 

Table  3  also  shows  the  standard  deviation  predicted  by  (22).  Note  the 
agreement  with  the  Monte-Carlo  simulations. 
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Table  4.  Sample  RMS  error  measurements,  1000  trials 


Method  Run 


Source  Location 


Source  Bearing  (xlfr) 


The  performance  of  the  SI  estimator  can  be  evaluated  by  comparison  to  the 
Cramer-Rao  (C-R)  lower  bound,  a  lower  bound  on  the  variance  of  any  unbiased 
estimator  [10].  Since  the  SI  method  is  biased,  the  root -mean-squared  (RMS) 
error,  defined  by 

a  -  /Bias2  +  Variance 

is  compared  to  the  C-R  bound.  Table  4  lists  the  Cramer-Rao  lower  bound 
standard  deviation  [<?]  along  with  the  RMS  error  for  the  SI  estimator.  As  seen 
in  this  table,  the  performance  of  the  SI  estimator  is  nearly  optimal, 
approaching  the  C-R  bound,  note  that  the  SI  estimator  appears  to  be 
relatively  unbiased,  as  the  RMS  error  is  nearly  equal  to  the  standard 
deviation  in  all  cases.  Therefore,  we  expect  (22)  to  accurately  predict  the 
RMS  error. 


The  extent  to  which  the  SI  method  is  not  an  optimal  least  squares  metnod 
is  a  function  of  the  extent  to  which  computed  by  (13)  is  not  equal  to 

;i x s a  computed  by  (14).  It  was  therefore  of  interest  to  compare  these 

quantities.  In  the  high-noise  case  (a,  =  1)  above,  for  both  source 
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positions,  the  quantity  h  .  a x, a /R  )  was  typically  less  than  O.nns  and 
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almost  never  greater  than  0.01.  Thus,  the  SI  method  is  very  close  to  an 
optimal  least-squares  equation-error  method  under  the  above  conditions.  It 
would  be  of  interest  to  find  out  when  and  if  the  disparity  between 

*■  A 

Rs  and  a_x_s a  becomes  large,  and,  if  so,  what  effect  this  has  on  the  source 
location  estimate. 


4.  Geometric  Interpretation 

This  section  gives  a  geometric  interpretation  of  the  SI  solution  and 
discusses  it  in  light  of  the  Monte-Carlo  results  presented  above. 

4.1  The  Error  Criterion 


"Is 


The  goal  is  to  localize  the  source  x  and  therefore  to  minimize 
—  — 

.  x  a  for  some  norm.  We  will  only  consider  l  (sum-of-squares )  norms. 

“5  2 


Since  only  the  RD  measurement  vector  d_  is  known,  a  natural  formulation  would 
appear  to  be 


Minimize  iuh  =  :]£  -  d_(x_ )  u  (23) 

x 

-  -  — s 

where  _d (x_  )  is  the  vector  of  all  measurable  RDs  corresponding  to  the  source 

^  A 

location  estimate  x_s  .  This  error  criterion  is  especially  well-suited  to  the 

case  of  zero-mean  errors  in  the  RD's.  Indeed,  if  the  RH  errors  are  Gaussian 
perturbations,  then  (23)  provides  the  maximum  likalinood  estimate  for 
x  (which  is  now  regaraed  as  a  parameter  determining  the  mean  d(x  )  of  the 
nultivariate  normal  distribution  for  d).  For  this  reason,  the  solution  to 
(23)  will  be  referred  to  as  the  maximum  likelihood  estimate  of  the  source 
1 ocation. 

As  mentioned  in  the  introduction,  solving  (23)  requires  nonlinear 
minimization  techniques.  For  this  reason,  the  SI  method  does  not  solve 
(23).  Instead,  it  approximately  minimizes  the  L?  .norm  of  an  "equation  error" 
which  was  chosen  purely  to  simplify  the  solution.  As  seen  in  the  Monte-Carlo 
results,  good  estimates  are  obtained  nonetheless. 

°lacing  sensor  I  at  the  origin  as  before,  and  using  only  N-i  3n's,  all 


15 


referred  to  sensor  1,  (23)  can  be  interpreted  as  finding  the  sphere,  passing 
through  sensor  1,  whose  surface  is  as  close  as  possible  to  being  d.^  away 
from  the  ith  sensor.  This  arrangement  is  shown  in  Figure  3  for  the  noiseless 
case. 

As  seen  in  Figure  a,  the  sphere  of  radius  R  ,  centered  at  x  and 

S  — 5 

passing  through  sensor  1,  is  a  surface  of  zero  RD  to  sensor  l.  The  distance 

A 

from  the  sphere  to  sensor  i  is  then  d^.  The  problem  is  then  to  position 
such  that  d.^  «  d  for  every  i.  Accordingly,  the  error  minimized  in  (23) 
is  equivalent  to  the  sum  of  squared  differences  between  the  sphere-to-sensor 
distance  d- ^  and  the  measured  RD  d.^  for  that  sensor. 

Since  the  sphere  around  the  source  must  always  pass  exactly  through 

sensor  i,  the  solution  is  sensitive  to  the  choice  of  sensor  l.  Improved 

results  in  the  case  of  a  systematic  bias  in  measuring  the  RD's  relative  to 

sensor  1  may  be  obtained  by  removing  the  constraint  that  the  source  sphere 

must  pass  through  sensor  1.  This  can  be  accomplished  by  adding  a  constant 

Y  to  the  RD  vector,  and  solving  for  y  as  well  as  R  and  x  . 

s  — s 

4.2  Spherical  Interpolation 


We  first  show  that  the  equation  error  (7)  (approximately  minimized  by  the 

SI  method)  is  closely  related  to  the  maximum  likelihood  error  minimized  in 

2 

(23).  Adding  and  subtracting  Rs  in  the  definition  of  the  equation  e^ror  (7) 
gives  (upon  introducing  hats  to  denote  estimated  quantities) 

2x[v  R*  -  (d^  +  2Rsdn  +  R*),  i  =  2 , 3 , . . . ,  N 


ei  =  :|*i  -  i1'2  -  <«s  *  du)2 


(24) 


where  x^  and  R?  are  the  estimated  source  location  and  range,  x.  is  the 

itn  sensor  location  (known  exactly),  and  d.  is  the  measured  ranae 

1 1 

difference.  Let 
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denote  the  RD  predicted  by  the  source  location  estimate  ^  (cf.  Figure  4). 
Then  (24)  can  be  written  as 

*t  ■  (»s  *atl)2.  (5S  *du)2 

=  »s<d11  '  d11>  -  dfi  -  afi  (25) 

■  <ail  *  dil  *  25s><du  -  dil> 

Assuming  that  the  noise  in  the  delay  estimates  is  small  compared  to  the  delay 
values,  e .  can  be  written  as 

•t  ■  2<»s  *dil>(dU  -  dH) 

This  form  of  the  error  displays  the  equation  error  as  the  maximum  likelihood 

A  A 

error  d.,  -  d.,  times  the  term  R  +  d...  When  the  source  range  is  large 
UU  s  1 1 

compared  to  the  intersensor  separation,  the  SI  equation  error  reduces  to 


ei  "  -Rs(di 1  -  dil} 

The  difference  between  the  SI  method  and  the  maximum  likelihood  method  is  the 
tendency  of  this  error  to  pull  the  estimated  source  location  toward  the 

A 

origin,  thereby  making  Rs  and  the  equation  error  smaller.  This  contraction 
of  the  source  estimate  toward  sensor  l  was  consistently  observed  in  the 
simulations.  Note  that  there  is  no  difference  between  the  maximum  likelihood 
solution  and  the  SI  method  with  respect  to  bearing  estimation  when 

R;  »  |dn|  <■  |Sn|. 

The  above  discussion  leads  naturally  to  a  weighting  function  and  an 
iterative  technique  for  obtaining  the  solution  to  (23).  The  weighting 
matrix,  W  ,  in  (II)  should  be  given  by, 


W  =  diagf- 


'dir<V2fV 


(26) 


where  diag(v.)  is  a  diagonal  matrix  with  v.  as  the  element  in  the  ith  row 
and  column,  and  q.,  +  d-,  2R$  is  given  by  a  prior  solution  of  (7),  or 

estimated  from  a  priori  information.  The  solution  of  (23)  is  then  found  by 


iteratively  solving  (7)  with  successively  updated  values  of  W  using  (26). 
Assuming  this  iteration  converges,  the  weighting  cancels  the  first  term  of 
(25)  and  the  maximum  likelihood  error  remains.  Conditions  for  the  convergence 
of  this  scheme  to  the  solution  of  (23)  need  to  be  determined. 

5.  Summary 

In  this  paper,  a  closed-form  solution  for  localizing  a  single  source  in 
n-dimensions  from  TDOA  information  gathered  from  an  N  >  n+l-sensor  array  was 
described  and  evaluated  by  Monte-Carlo  simulations.  It  was  found  that  the  SI 
method  exhibited  an  RMS  localization  error  close  to  the  Cramer-Rao  lower 
bound.  Finally,  in  support  of  our  simulation  results,  we  showed  that  the  SI 
method  is  closely  related  to  the  maximum  likelihood  estimate  for  the  case  of 
Gaussian  TDOA  measurement  errors. 
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Figure  Captions 


Figure  1.  Diagrams  illustrating  notation  and  certain  geometric  relations  for 

the  case  of  a  single  source  x  .  Labels  imbedded  within  a  line 

*— $ 

denote  the  length  of  the  corresponding  vector.  For  example, 


Figure  2.  Source-sensor  geometry  used  in  all  simulations. 


Figure  3.  Sensor  array,  source  and  Rn  geometry.  The  sphere  of  radius 

R$  drawn  around  the  source  is  the  surface  of  zero  RD  relative  to 
sensor  1.  The  perpendicular  distance  from  the  sphere  to  any  sensor 
is  the  RD  for  that  sensor  relative  to  sensor  1. 


Figure  4.  Geometric  representation  of  the  relationship  given  in  equation  (7). 
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Abstract.  Cramer-Rao  lower  bounds  are  derived  for  the  variance  of  unbiased  estimators  of  source  location, 
source  range,  and  source  direction  given  range-difference  (pd)  measurements.  The  Cramer-Rao  bound  (orb) 
for  range  estimate  variance  is  found  to  be  a  quartic  function  of  source  range.  The  '’PC  for  bearing  estimate 
variance  is  found  to  be  insensitive  to  the  source  location,  an  inverse  quadratic  function  of  the  sensor-arrav 
solid-angle  as  seen  from  the  source,  and  an  inverse  linear  function  of  the  spatial  density  of  the  sensors  within 
the  array  (assuming  independent  rd  measurements  in  each  pair  of  sensors).  The  theoretical  bounds  are 
compared  to  empirical  mean-squared  errors  obtained  using  a  fast,  suboptimal.  least  squares  estimator  of 
source  location  called  the  Spherical  Interpolation  (si)  method  [1|.  The  si  method  estimates’  mean-squared 
errors  are  found  to  be  in  the  range  of  1.0  to  1.5  times  the  respective  Cramer-Rao  lower  bounds. 

1.  Introduction 

The  problem  of  locating  an  object  by  analyzing  signals  received  from  it  is  a  basic  one  in  the  fields  of 
underwater  acoustics,  geophysics,  aerospace  and  industry.  It  is  therefore  of  interest  to  know  how  well  a 
source  localization  algorithm  can  perform  when  the  received  signals  are  corrupted  by  noise. 

A  source  tells  of  its  location  through  the  relative  time  delays  measured  among  the  signals  it  radiates  to 
the  sensors  of  an  array.  In  a  constant-velocity  medium,  the  time-dilference-of~  arrival  (tD"a)  between  signals 
received  at  two  sensors  is  proportional  to  the  difference  in  source-sensor  range,  termed  range  difference  (r? D ), 
and  places  the  source  on  a  hyperboloid  of  revolution  about  a  line  drawn  between  the  sensors.  In  p  dimensions. 
HD  measurements  from  p  non-coplanar  sensor  pairs  [p  1  sensors)  are  sufficient  to  localize  a  source  l’2l. 

Typically,  the  no’s  are  not  known  precisely  and  the  source  location  has  to  be  estimated.  The  source 
location  estimate  will  have,  under  suitable  assumptions,  an  associated  variance  and  bias.  The  Cramer-Rao 
bound  (•  no)  gives  a  lower  bound  on  the  variance  of  an  unbiased  estimator,  and  is  often  used  in  evaluating 
the  performance  of  an  estimator.  An  unbiased  estimator  which  achieves  the  rd  is  called  efficient;  it  is  well 
known  that  all  Maximum  Likelihood  (ml)  estimators  are  efficient.  Unfortunately,  in  the  case  of  estimating 
source  location  from  noisy  pd  s,  the  ml  estimator  is  expensive  to  implement,  requiring  nonlinear  optimization 
techniques.  While  much  is  known  about  the  problem  of  estimating  Pd’s  from  received  signals  (a  recent  special 
issue  on  the  topic  is  l13l),  there  has  been  relatively  little  work  on  the  estimation  of  source  location  from  RD 
measurements  j  1 , 2 . 3 , 4 , 5 j . 

This  paper  presents  Cramer-Rao  lower  bounds  on  the  variance  of  unbiased  estimates  of  source  location, 
direction,  and  range  from  RD  measurements.  These  bounds  are  then  compared  to  the  measured  variances 
of  estimates  obtained  using  the  Spherical  Interpolation  (si)  method  l!.  The  -i  method  provides  a  fast, 
closed-form,  least-squares  type  solution  for  source  location  from  RD  measurements,  and  is  shown  here  to 
have  rms  error  approaching  t. he  Cramer-Rao  lower  bound. 

The  organization  is  as  follows.  In  section  2  the  Spherical-Interpolation  Method  of  source  location 
estimation  is  derived.  In  section  3  the  Cramer-Rao  lower  bounds  on  the  variance  of  source  location,  bearing, 
and  range  estimates  are  derived  and  discussed.  Next,  section  4  presents  Monte-Carlo  simulation  results 
demonstrating  the  near-efficiency  of  the  -’I  technique. 

2.  The  Spherical- Interpolation  Method  for  Closed-Form  Localization 

Let  .V  denote  the  number  of  sensors,  and  let  i,,  denote  the  rd  between  sensors  r  and  ;  (;,  '  -  1 . .V). 

The  vector  of  (z.  y.  r)  spatial  coordinates  for  the  ith  sensor  is  denoted  z, ,  and  the  position  of  the  source  is 


i„j  iji  i].  i.l  i  j.  v.'-"--'-  ■“-  *. 


denoted  £,.  The  distance  between  the  source  and  sensor  i  is  denoted  by  D,  =  l|i,  -  i.  i|,  and  the  distance 
from  the  origin  to  the  point  x,  is  denoted  Rj.  Similarly,  f?.  =  l|x„  i|.  These  quantities  appear  in  Figure  la. 
Note  that 

d„  =  £,  -  D},  i  =  1, . . . ,  ;V,  j  =  l . N  (1) 

The  localization  problem  is  to  determine  x,  given  d ,,  for  i  and  ;  between  1  and  ;V  Note  that  there  are 
N{N  —  l)/2  distinct  rd's  dt,  (excluding  t  =  and  counting  each  d.,  =  -d;,  pair  once);  however,  any  W  —  1 
RD  measurements  which  form  a  "‘minimal  spanning  subtree’  determine  all  the  rest  (in  the  noiseless  case). 
The  redundancy  of  the  complete  set  of  rd  measurements  is  used  to  increase  noise  immunity. 

2.1.  Equation- Error  Formulation 

To  solve  the  localization  problem,  we  first  map  the  spatial  origin  to  an  arbitrary  sensor,  say  the  ;th, 
i.e.,  xt  —  0.  The  rd  d,,  then  reduces  to 


d,,  =  il  x,  —  x.  !j  -  '|  x,  '| 

and,  from  the  definitions  above, 

d„  =(/?M?--2xJi,)i  -R.  (2) 

Adding  R.  to  both  sides  of  (2)  and  squaring  gives 

(d, J  +  R.)2  =  R:  -R; '  -  2x^x.  (3) 

Equation  (3)  is  the  Pythagorean  theorem,  illustrated  in  Figure  lb.  Moving  all  terms  to  the  right  side  of  (3), 
the  /?;  terms  cancel,  and  we  are  left  with 

0  =  R;-l;,  -  1RJ„  -  2x;x.  (4) 

The  ;th  equation  is  degenerate  so  we  have  iV  -  1  equations  in  three  unknowns  x,. 

As  the  delays  are  not  known  precisely,  we  introduce  a  so-called  “linear  equation  error’  81  into  the 
right- hand-side  of  (4),  and  minimize  it  in  a  least  squares  sense  to  provide  an  estimate  of  the  true  solution. 
Without  loss  of  generality,  let  ;  =  1.  Then  (4)  becomes 

*.  =  R?  ~  dfi  ~  -d?.d, i  -  2x,Tx. ,  t  a  2.  3 . ,V  (5) 


where  c,  is  the  linear  equation  error  to  be  minimized.  The  set  of  .-V  -  1  equations  (5)  can  be  written  in 
matrix  notation  as 


e=o-  IR.d  -  2Sx 


(«) 


where 


R2-.at  , 

/  dq  i  \ 

/  Uq  xq  \ 

Rl-Ai 

,  A 

d,i  1 

xi  Va  -a 

; 

,  4  = 

•  1 
• 

S  = 

'  Ry  ~  d  v  i  ^ 

d:Vl  1 

'  Jt.v  y.v  x.v  ' 

It  is  worth  noting  that  equation  f <>)  is  line ar  in  x,  given  R.,  and  it  is  linenr  in  R,  given  r.  Error  vectors 
which  are  linear  in  the  unknowns  yield  closed-form  least  squares  solutions. 

The  least  squares  solution  for  x.  given  R .  is 


X,  =  ,s;vl>  -  2  f? .  i) 

where  the  weighted  Iine.tr  equation  error  energy  Jfx.j  =  -:W-  is  minimized  for 

S{v  =  (STWS)'‘STW 


(3) 
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and  W  is  symmetric  positive  definite  (or  simply  diagonal  and  positive).  If  W  =  I,  S,Y  becomes  the  pseudo¬ 
inverse  of  S. 

To  obtain  a  true  least-squares  estimate,  it  is  necessary  to  minimize  -/(?.)  with  respect  to  x.  while 
allowing  R.  to  vary,  maintaining  the  relation  R .  =  ![z,  j|.  This,  unfortunately,  is  a  nonlinear  minimization 
problem,  and  serious  difficulties  can  occur  when  J(x.)  is  not  a  convex  functional  of  x,.  However,  the 
nonlinearity  can  be  eliminated  as  described  in  the  next  section. 


2.2.  The  Spherical  Interpolation  Method 


The  basic  idea  of  the  si  solution  is  to  substitute  (8)  into  (6)  and  minimize  the  linear  equation  error 
again,  this  time  with  respect  to  R..  This,  surprisingly,  yields  a  linear  least  squares  problem  for  finding  f?,, 
and  the  solution  is  quick  and  inexpensive.  The  technique  is  made  possible  by  the  fact  that  the  formal  least 
squares  estimate  of  x.  given  R.  in  (6)  is  itself  linear  in  R..  When  the  minimizing  R.  value  is  found  in  this 
new  linear  equation,  the  corresponding  value  of  r,  (via  (8))  is  automatically  a  minimizer  of  the  squared 
equation-error  norm  with  respect  to  x,  given  this  R.. 

Rewriting  the  linear  equation  error  defined  in  (6)  to  eliminate  x„  by  subsituting  the  value  from  (8)  yields 


-6  -  2R.d  -  SS,V(o  -  2R,d)  =  (I  -  SS^)(d  -  2 R.d) 


where  I  is  the  .V  —  1  by  ,V  —  l  identity  matrix.  Now,  *  is  linear  in  the  single  unknown  R..  Define  the  ,V  -  1 
by  iV  -  1  matrices 

P>  =  SSiV  =  S(STWS)-‘STW 


(10) 


Pt  =  I  -  SS,V  -  I  -  P- 

The  rank  of  P.  is  at  most  3  regardless  of  its  order  ,V  -  i,  and  Pi-  has  rank  at  least  (jV-  1)  -3.  Also,  P-  and 
Pt  are  idempotenc,  i.e..  P*  =  P-  and  (Pi)'  =  Pt  Idempotent  matrices  can  be  interpreted  as  projection 
operators.  For  example,  when  W  =  I,  the  operator  P,-  projects  an  ( i V  —  l)-vector  into  the  subspace  spanned 
by  the  columns  of  S,  and  the  operator  Pi  projects  into  the  orthogonal  complement  of  the  subspace  spanned 
by  the  columns  of  S.  Thus,  for  example,  P.x  is  orthogonal  to  Pty  for  every  pair  of  vectors  x  and  y. 

In  the  uondegenerate,  four-sensor  case,  Py  =  /3,  and  the  error  £  is  zero.  In  the  more  general  case  of  .V 
sensors, 

1  -  Pt  (£  —  -R.d) 


so  that 


|  e  ilfr  =  iTW£  =  (d  -  2/?,d)TPiWPi(£-2f?,<i)  =  (£  -  2R,d)T  W,x  (d  -  2R.d) 


Minimizing  J'  with  respect  to  R,  is  a  form  of  weighted  least  squares  in  which  the  weighting  matrix  W,i  is 
of  rank  N  —  4.  The  missing  dimensions  reflect  the  degrees  of  freedom  removed  by  choosing  sensor  1  as  the 
origin  and  substituting  in  the  least-squares  solution  (8)  for  the  three  spatial  source  coordinates.  Note  that 
when  W  =  I.  =  P*.  The  least-squares  minimizer  of  J'{R.)  is  given  by 


2  irW..d 

Substituting  this  solution  into  (8)  yields  the  closed-form  source  location  estimate 


(ID 


d) 


(12) 


Clearly,  the  computational  burden  of  (12)  is  very  low  compared  to  iterative  nonlinear  minimization. 

Note  that,  in  general.  R,  as  computed  by  (11)  is  not  necessarily  equal  to  1  '  computed  from  (12) 

except  as  the  ftp  noise  approaches  zero.  Therefore,  we  define  the  range  estimate  by 


/?.,  4  :|x,, 


[13! 
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Similarly,  the  bearing'  estimate  is  defined  as  the  vector  of  direction  cosines  from  the  origin  at  sensor  1  to 
the  source: 

&/*«#*!  (“) 


3.  The  Cramer-Rao  Lower  Bound 

The  Cramer-Rao  bound  (« jr B )  is  a  lower  bound  on  the  variance  of  an  unbiased  estimator  *9, 10, 11.121; 
estimators  achieving  this  bound  are  said  to  be  efficient.  In  this  section  •  ’n B ’5  are  derived  for  unbiased 
estimators  of  source  location,  source  range,  and  source  direction  based  on  rd  information. 

We  first  review  the  ''RD,  the  reparameterization  relation,  and  the  information  inequality.  The  variance 
of  estimating  the  tth  element  of  9_  with  an  unbiased  estimator  9,  based  on  an  observation  vector  z  is  bounded 
below  by  the  ith  diagonal  element  of  the  inverse  of  the  Fisher  information  matrix  .91: 


Var(0.)  =  E  {*7}  -  E2  {(?.}  >  [a,-1]  _ 


where  Var(  )  denotes  variance,  E{  }  denotes  expectation,  (M|tl  denotes  the  tth  diagonal  element  of  the 
matrix  M,  and  is  the  Fisher  information  matrix  for  9,  defined  by  [9| 


[f  9  is  parametrized  by  t  (i.e.,  9  =  hit)),  the  Fisher  information  matrix  for  t  can  be  written  in  terms  of  the 
Fisher  information  matrix  for  9_  via  the  reparameterization  relation  ;9| 

-  I  -  Ol\ 
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where,  the  (t,  ;)th  element  of  d9J dt  is  d9;/dtr  Finally,  if  7  is  an  unbiased  estimator  of  the  scalar  •?(£).  the 
variance  of  7  is  bounded  by  the  information  inequality: 


?lY  dJ. 


Var(q)  > 


We  now  derive  the  •  Rb’s  for  source  location,  range  and  direction.  Defining  z.(t)  as  the  tth  element  of 
z„,  the  Cramer-Rao  lower  bound  on  the  variance  of  an  unbiased  estimator  of  source  location  z„,  based  on 
RD  estimates  d,  is  given  by  (15)  as 

Var(x.(t))  >  (18) 

where.  ~T  is  the  Fisher  information  matrix  for  z, .  Using  the  information  inequality  (17),  the  variances  in 
estimating  source  range  and  source  direction  from  rd  measurements  are  bounded  below  by 


where  n,(i)  is  the  ith  element  of  Q,. 

Since  the  source  location  r,  is  a  parameter  in  determining  the  distribution  of  the  RD  measurement 
vector  i.  by  the  reparameterizatioii  relation  (L6),  :r  can  be  written  as 


4 


and.  as  shown  in  Appendix  A.  assuming  the  RD  estimates  are  Gaussian  distributed  with  correlation  matrix 
R.i,  assumed  independent  of  source  location,  F.i  =  R//1  and 


It  should  be  noted  that  when  the  Rd's  are  estimated  from  Gaussian  signals,  commonly  used  estimators  are 
unbiased  and  have  j6,7| 

R.i  as  <y\l  (22) 

where  <y\  is  the  variance  of  each  individual  P.D  and  is  dependent  only  on  signal  and  medium  characteristics. 

To  find  the  Cramer-Rao  lower  bound  on  the  variance  of  the  source  location  estimate,  it  remains  only 
to  find  the  sensitivity  (derivative)  of  range  difference  with  respect  to  source  location  dd/dx..  Similarly,  the 
RB  s  for  range  and  bearing  of  the  source  are  easily  calculated  after  finding  the  sensitivities  of  range  and 
bearing  with  respect  to  source  location. 

3.1.  Sensitivity  of  Range  Difference  to  Source  Location 

The  derivative  dd/dx.  can  be  evaluated  by  differentiating  equation  (6)  in  the  noise-free  case  with  respect 
to  x.  and  solving  for  dd/dx,.  Setting  s  =  0  in  (6)  and  differentiating  with  respect  to  x,  yields 

dS  d{R.d) 

—  2  —7 — —  —  -S  =  0  (23) 

dx,  dx.  dx. 

From  (7),  the  ith  element  of  S  is  R~  —  dfL  so  that  dS_dx,  becomes 


where.  A i  is  a  diagonal  matrix  with  d,  =  <i, t  as  the  tth  row/column  entry.  The  derivative,  d{R.d)/dx,,  is 
evaluated  using  the  product  rule: 

d[R.d)  3d  3R.  3d 

-  =  7?.  t - - iQ.  (25) 

dx.,  dx,  dx,  dx, 

where.  R.  =  ‘|r.  !  has  been  used  to  evaluate  3R.;dx.  =  Q,  =  r./'lx.  |.  Substituting  (24)  and  (25)  into 
(23)  with  dx..  dx,  -  I  yields  dd/dx.: 

A^^A=0  ~  p-  =  -VA  (26) 

dx,  dx, 

where 

V  =  A-1 

A  =  Aj-  R.I  (27) 

A  =  dQj  *  S 

Post-multiplying  through  (6)  by  Q J//?.,  solving  for  d  Q] .  and  substituting  into  (27)  produces  a  formula 
for  A  in  terms  of  £  rather  than  d: 

'V.T 

A  =  SP7-,fl'P,  (23) 


p£_  =  x.fxjx.r^r  =  n.n: 

P;  =  (I  -  Pr.) 

The  idempotent  matricies  P-  and  Pc  are  projection  operators  which  remove  the  vector  components  per¬ 
pendicular  to  x.  and  parallel  to  r,,  respectively.  Thus,  A  multiplies  the  :d  subspace  orthogonal  to  Q.  by 
the  matrix  S.  and  it  ■‘bends'’  the  id  component  parallel  ro  Q,  into  the  direction  of  \  scaled  by  j  ■'  1/27?.. 


5.2.  Cramer-Rao  Bounds  for  Location,  Range,  and  Bearing 

The  Cramer-Rao  lower  bound  on  the  variance  of  an  estimator  of  source  location,  x,,  from  Gaussian 
distributed  RD  measurements  is  obtained  by  using  (26)  in  (21): 


Var(x,  (i))  > 


=  (atv/?.71va) 


»cr.i(ATV3A)" 

The  Cramer-Rao  lower  bound  on  the  variance  in  estimating  source  range.  /?.,  from  RD  measurements 
is  similarly 

VI*.)  *  (£)>-(£)’.  or/ro. 

=  q: (ATvf?;‘vA)  ‘a,  (31) 

*  cr-tQ!  (ATV:A)~‘  Q. 

For  the  bearing  estimate  'RD,  we  need  (using  Q.  =  z,,  ']  x,  1): 


1  P' 

_  it  n  -  -■ 


=  i-~—  =  — (i  -  n.ni)  =  ~ 

(xJx.)J  R.  R- 


The  ro  for  the  bearing  escmiate,  Q,,  is  thus 


dx.)'r-  \9x.)  R-  -•  -■ 

=  (a7vr;1va)'‘p; 

*  ^P;  (A7v: A)”1  P; 


5.5.  Far- Field  Analysis 

In  this  section,  we  analyze  the  location,  range,  and  bearing  pc’s  in  the  “far  field'’  (/?.  i.^ri) 
and  the  '‘infinite-range  case”  (R.  —  oo).  These  assumptions  on  the  scale  of  R.  greatly  simplify  the  rb 
expressions  and  give  insight  into  the  behavior  of  the  bounds.  We  only  consider  the  i.i.d.  case  R.i  -  c-I.  as 
is  often  done  in  practice  6.7j. 

The  far  fold 

The  far  field  of  the  sensor-array  aperture  is  defined  by 


R.  »  d.  i . 


.V 


At  this  range,  we  have  A  =  R, I  —  A.j  ss  R.l.  and  from  (22),  (27),  and  (20) 


7-1 

•  T. 


7^:(AtA)"1  (36) 

Using  (23)  for  A  and  (31)  for  the  -no  variance  of  R,,  (36)  gives  die  far-fieid  range-variance  bound  as 
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Var(  R. )  ~  4a]  R4. 


L£Tptj£ 


where 


Pu  =  i-  ufirur'u1 

U  =  SO; 

nt  = ;  Qr  Q?  ’ 


(38) 


where  Qj"  and  07  are  any  two  mutually  orthogonal  vectors  in  3-space  which  are  also  orthogonal  to  Q_.  The 
matrix  P?r  is  the  idempotent  projection  into  the  orthogonal  complement  of  the  subspace  spanned  by  the 
columns  or  U  =  Sflt 

Similarly,  the  far-field  bearing-variance  bound  is  approximated  by 

Var(fi,(t))  ~  a]nt  (uTp,_u)~1  ntT| 


_  o  o 

P,*  =  1  -  =5- 

6  ‘  o 


(39) 


and  again  U  =  Sflt 

We  now  discuss  qualitative  relationships  revealed  by  the  above  formulas.  Equations  (37)  and  (39)  show 
that  the  far-field  range  estimate  Cramer-Rao  lower  bound  variance  increases  quartically  with  the  range,  and 
the  far-field  bearing  escimate  deviation  is  constant  with  source  range.  Consequently,  the  minimum  standard 
deviation  of  a  source  location  estimator  increases  with  the  square  of  the  source  range  in  the  source  direction, 
and  increases  linearlv  with  the  range  in  the  other  directions. 

Second  from  the  far-field  bounds,  it  is  apparent  that  increasing  the  norm  of  i\  reduces  the  Cramer- 
Rao  bound  variance  in  estimating  source  location,  direction,  and  range.  By  the  definition  of  A,  given  in 
(28),  therefore,  the  far-field  bounds  are  reduced  when  the  size  of  S  in  the  directions  perpendicular  to  the 
source  direction  increases,  or  when  the  norm  of  i  increases.  The  sizes  of  S  and  6  are  increased  by  placing 
the  sensors  far  from  each  other  and  far  from  sensor  1,  as  seen  by  the  source,  i.e.  by  increasing  the  size  of  the 
array.  The  norm  of  e  is  also  increased  by  aligning  the  array  such  that  lines  drawn  between  sensor  1  and  the 
other  sensors  tend  to  be  perpendicular  to  the  source,  thereby  decreasing  d  and  increasing  o.  Increasing  S 
or  '  lowers  the  source  estimate  variance  by  increasing  the  effective  aperture  size  of  the  sensor  array,  as  seen 
from  the  source. 

Third,  the  variances  in  the  source  location,  direction,  and  range  estimates  are  all  proportional  to  the 
?D  noise  variance. 

To  corroborate  the  above  observations,  the  exact  Cramer-Rao  bounds  were  evaluated  numerically  and 
compared  to  the  far-field  observations.  Figure  2a  shows  the  square  root  of  the  ’RD  range  estimate  standard 
deviation  plotted  against  range  for  the  case  of  a  100-iueter-wide,  200-meter-long  9-sensor  prism-shaped  array 
locating  a  source  placed  at  various  ranges  along  a  line  through  the  origin  (see  Monte-Carlo  Run  1  in  Table 
1).  The  range  estimate  standard  deviation  can  be  seen  to  increase  very  linearly  with  the  square  of  the  range, 
even  at  small  ranges.  Figure  2b  shows  the  rd  direction-cosines  estimate  standard  deviation  plotted  against 
ranee  for  the  case  described  above.  The  bearing  estimate  standard  deviation  can  be  seen  ro  he  constant  with 
increasing  range,  except,  ac  small  ranges  where  there  appears  to  be  a  point  of  maximum  angular  resolution, 
or  best  focus. 
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The  Infinite- Range  Case 


In  the  infinite-range  case.  R. 
(6)  gives 


—  oo,  only  the  source  direction,  Q.  is  of  interest. 

-sn,  =  d  =>  n,  =  -s‘i 


At  infinite-range  equation 
(40) 


where  S"  =  (STS)  lST  Note  that  by  this  defining  relationsip,  the  infinite-range  Q,  is  not  constrained  to 
have  unit  norm. 

The  Fisher  information  matrix  for  the  infinite-range  Q,  estimate  is  easily  computed  from  *4  =  R  t  1  and 
the  reparameterization  relation;  accordingly  the  infinite-range  Cramer-Rao  bound  for  the  henring  estimate 


is 


where 


=  S'TZ?.,S* 


*^(STS) 


(411- 


Note  that  this  bound  is  independent  of  source  direction  (since  !|  0.  j|  is  unconsrtained).  It  is  also  inversely 
proportional  to  the  square  of  the  scale  of  the  array,  i.e.,  if  the  array  were  twice  as  large,  the  ■  imj  variance  in 
estimating  Q,  would  be  reduced  by  a  factor  of  four.  The  bound  is  also  inversely  proportional  to  sensor-array 
density;  e.g.,  if  the  number  of  sensors  is  doubled  by  measuring  two  independent  HD  estimates  at  each  sensor 
(thus  doubling  the  length  of  S  by  repeating  each  row),  the  value  of  STS  doubles  and  the  variance  of  r, 
halves.  In  practice,  the  bound  is  inversely  proportional  to  array  density  up  to  the  density  limit  for  which 
the  no’s  can  be  assumed  uncorrelated.  As  always,  the  ■  RD  variance  is  proportional  to  the  no  variance, 
it  is  shown  in  Appendix  B  that  the  infinite-range  lower  bound  (41)  is  achieved  by  the  >t  method. 
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4.  Monte-Carlo  Simulation  Results 

This  section  presents  Monte-Carlo  simulation  resuits  comparing  the  performance  of  the  -I  method  to 
tht  Cranier-Rao  lower  bound.  The  simulations  were  implemented  in  the  Ctrl-C*  language  on  a  Vax  11  735 
computer.  Since  the  estimator  exhibits  a  slight  bias,  the  root-inean-squared  (VimsJ  error,  defined  by, 

i7 r  =  \E(x.  -  x.):l  - 

of  the  -I  method  was  compared  to  the  Cramer-Rao  lower  bound  standard  deviation.  It  was  found  that  the 
-t  method  is  nearly  efficient  in  all  cases  tried. 

The  simulations  employed  the  two  nine-sensor  arrays  shown  in  Figure  3,  with  the  two  source  locations 
shown  and  two  choices  of  sensor  1.  The  i?D  measurements  were  simulated  by  adding  i.i.d.  white  Gaussian 
pseudo-random  noise  to  the  true  rd  values  (corresponding  to  the  practical  assumption  R.±  =  o’" I  1.6,7'). 
Sample  bias,  variance,  and  RMS  error  of  the  unweighted  .~r  method  were  computed  by  averaging  the  results 
of  100-trial  Monte-Carlo  runs.  The  Cramer-Rao  bounds  for  each  Monte-Carlo  run  were  computed  using  the 
formulas  given  in  section  3.2  above. 

Table  1  describes  the  environmental  information  (source  location,  range,  bearing,  additive  r?D  noise 
level,  and  sensor  array)  for  each  of  the  eight  100-tna!  Monte-Carlo  runs. 

Tables  2.  3  and  4  show  the  sample  pms  error  of  the  m  source  location,  range  and  bearing  estimates 
and  the  corresponding  Cramer-Rao  bounds  for  each  of  the  Monte-Carlo  runs.  It  should  be  noted  that  the 
sample  biases  of  the  ,-i  estimates,  not  shown  here,  were  typically  less  than  one  tenth  of  the  sample  standard 
deviation,  and  were  never  prominent  in  determining  the  P  M-  error. 

From  Tables  2.  3,  and  4.  as  well  as  other  results  not  presented  here,  the  -I  -stimarors  of  range,  bearing, 
and  source  location  are  seen  to  have  pms  error  in  the  range  1.0  to  1.5  times  the  amer-Rao  lower  bound 
standard  deviation.  The  -[  estiniator  is  therefore  performing  almost  efficiently,  in  other  terms,  the  -t 
estimator  approximates  the  maximum  likelihood  estimator  of  source  location  given  no  measurements. 


3 


J 


*  Ctrl-C  is  a  Mat  Lab- based,  high-level,  ‘‘matrix  calculator’  language  sold  by  SCT,  Inc. 
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It  is  interesting  to  see  that  the  choice  of  sensor  1  (the  sensor  to  which  the  fd’s  are  measured)  has  a 
large  effect  on  the  Cramer-Rao  lower  bound  (as  much  as  a  factor  of  two  in  standard  deviation).  It  has 
been  observed  that  the  range,  and  therefore  source  location  estimates,  have  smaller  ■  tb  for  sensor  1  choices 
which  are  either  centrally  located  in  the  sensor  array  relative  to  the  source  or  are  closer  to  the  source.  Also 
observed  is  the  tendency  for  the  '"F  b  to  be  smaller  for  source  direction  estimates  for  which  sensor  1  is  chosen 
to  be  centrally  located  in  the  sensor  array  as  seen  from  the  source. 

The  above  observations  are  consistent  with  the  far-field  analysis  of  section  3.  For  these  sensor  1  choices, 
R,  is  smaller  and  ||  0  !|  is  larger,  both  resulting  in  smaller  Cramer-Rao  lower  bounds. 

Also  consistent  with  prior  analysis  is  that  the  ''re’s  obtained  using  the  larger  (prism  shaped)  array,  are 
lower  than  the  corresponding  bounds  using  the  smaller  (cube  shaped)  array.  By  the  analysis  of  Section  3, 
when  the  array  size  increases,  the  matrix  S  is  larger,  and  the  orb  standard  deviation  is  decreased. 

The  .-I  technique  is  seen  to  achieve  the  Cramer-Rao  lower  bound  for  some  choices  of  sensor  1,  while 
performing  at  as  much  as  1.5  times  the  ''FB  for  other  choices  of  sensor  1  for  the  same  source  location,  FD 
noise,  and  sensor  array.  It  has  been  observed,  however,  in  all  cases  tried,  that  the  choice  of  sensor  1  giving 
the  lowest  ''FB  also  gives  the  lowest  RMS  error  estimates  using  the  si  technique.  This  suggests  a  three  step 
procedure  for  locating  a  source:  first,  use  an  arbitrary  sensor  1  to  locate  the  source  using  the  <1  method; 
second,  using  this  source  location  estimate,  find  the  best  choice  for  sensor  1  using  the  Cramer-Rao  lower 
bound  formulas  above;  and  third,  use  the  st  method  to  locate  the  source  using  RD  measurements  made 
from  the  best  sensor  1.  A  final  refinement  would  be  to  perform  Gauss-Newton  descent  on  !jd  -  d(x.)  (|*  to 
converge  to  the  optimal  solution. 

5.  Summary 

In  this  paper  we  presented  Cramer-Rao  bounds  on  the  minimum  variance  in  estimating  source  location, 
direction,  and  range  from  fd  measurements.  It  was  found  that  the  •  rb  range  estimate  variance  was  a  quartic 
function  of  the  source  range,  and  sensitive  to  the  sensor  array,  and  the  choice  of  reference  sensor.  The  source 
direction-cosines  estimate  variance  was  seen  to  be  sensitive  to  the  sensor  array  used,  but  insensitive  to  the 
source  location. 

The  -1  estimator  of  source  location,  direction,  and  range  was  described,  and  was  shown  to  be  asymp¬ 
totically  efficient  (as  the  range  approaches  infinity)  in  the  case  of  Gaussian  distributed  range-difference 
measurement.  Additionally,  the  si  method  was  shown  by  computer  simulation  to  a  be  in  the  range  of  1.0  to 
1.5  times  the  Cramer-Rao  lower  bound  for  a  several  near-field  source  locations. 


Appendix  A  —  Source  Location  Fisher  Information  Matrix  for  Gaussian  Distributed  Range- 
Difference  Measurements 


In  this  appendix,  we  derive  the  relationship  (21). 

Let  p(d|  r.)  denote  the  probability  density  function  of  i  given  z.,  where  d  is  a  random  ( :V  -  l)-vector, 
and  z,  is  a  vector  of  dimension  3  which  parametrizes  the  distribution. 

In  the  Gaussian  case  we  have  ( 14) 


P(d j  X.) 


1  - 


where  i  =  E{<j)  is  the  mean  (a  function  of  x.)  and  R.t  =  E{djT}  -  idr  is  the  covariance  matrix  of  the 
random  vector  i  (assumed  independent  of  x,). 

The  Fisher  information  matrix  is  defined  by  9j 
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The  Cramer-Rao  bound  for  an  unbiased  estimator  of  the  parameter  x.  based  on  d  observations  is 


Var(rJ  >  »“l  = 


'3x,  /  -  \<9x. 


Appendix  B  —  Efficiency  of  the  Infinite-Range  SI  Bearing  Estimator 

For  an  estimator  to  be  efficient,  it  must  be  unbiased  and  have  variance  equal  to  the  Cramer-Rao  lower 
bound.  Here,  we  show  that  the  unweighted  si  estimator  of  Q„  is  efficient  when  die  source  is  infinitely  far 
away  from  the  sensor  array,  assuming  it  is  based  on  unbiased  pd  estimates  d  with  variance  R.i.  (Efficiency  of 
the  bearing  estimator  is  defined  relative  to  the  RD  distribution.  Therefore,  the  underlying  RD  estimator  itself 
need  not  be  efficient.  It  is  also  possible  to  allow  the  RD  estimator  to  be  biased  as  long  as  it  is  understood 
that  in  this  case  an  unbiased  bearing  estimate  is  one  whose  mean  is  the  bearing  corresponding  to  the  mean 
HD  vector;  in  this  case,  the  functional  dependence  of  the  biased  RD  mean  on  source  bearing  must  be  d  .-fined 
and  differentiable.  VVe  prefer  to  assume  an  unbiased  RD  estimator  since  the  added  complexity  of  allowing 
RD  bias  does  not  appear  to  add  significant  practical  value.) 

We  first  show  that  the  si  estimate  Qs[  of  the  infinite-range  source  direction-cosines  vector.  From  (12). 
in  the  limit  as  R.  — *  ao. 


n„7  =  ‘i-  =  -S'd, 

~  !  R. 


R.  —  oo 


where.  S'  =  (STS)  ‘ST  Taking  expected  values, 


E(Qsr)  =  Ei-S'd)  =  -S'E(d)  =  -S'd 

by  the  linearity  of  expectation  and  the  fact  that  the  RD  estimate  is  unbiased. 
The  variance  of  the  ~I  estimator  is  given  by 

var(Qs/(t))  =  [r(6.v, Qlr)  -  £(6>-/)£(iW)T]  r 

Substituting  for  Q.-i<  we  obtain 

Var(Q.w/(0)  =  [S*  ffiT(diT)  -  Efi)  £(d)Tj  S'T] 

wliich  is  the  Cramer-Rao  lower  bound  given  in  (41). 
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Figure  Captions. 


Figure  1.  Diagram  illustrating  notation  and  certain  geometric  relations  for  the  case  of  a  single  source  x.. 
Labels  imbedded  within  a  line  denote  the  length  of  the  corresponding  vector.  For  example,  D,  »  l|  x,  -  x,  j|. 

Figure  2a.  The  square  rooc  of  the  Cramer-Rao  bound  standard  deviation  for  estimating  source  range 
plotted  as  a  function  of  source  range  for  the  source  direction,  sensor  array,  and  rd  noise  level  given  in  Run 
1.  Table  1. 

Figure  2b.  The  Cramer-Rao  bound  standard  deviation  for  estimating  source  direction  cosines  plotted 
as  a  function  of  source  range  for  the  source  direction,  sensor  array,  and  no  noise  level  used  in  generating 
Figure  2a. 

Figure  3.  Source-sensor  geometry  used  in  all  simulations. 
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Table  1 

Monte-Carlo  Runs 


Run 

Number 

# 

Sensor  1 

(x,y,z) 

meters 

Source 

Location 

(x,y,z) 

meters 

Array 

# 

1 

(0,0,0) 

(160,390,170) 

1 

2 

(0,0,100) 

(160,390,170) 

1 

3 

(0,0,0) 

(160,390,170) 

2 

4 

(50,50,50) 

(160,390,170) 

2 

5 

(0,0,0) 

(820,350,560) 

1 

6 

(0,0,100) 

(820,350,560) 

1 

7 

(0,0,0) 

(820,350,560) 

2 

8 

(50,50,50) 

(820,350,560) 

2 

Sensor  Array  1:  0  0  0 

0  0  100 
0  0  200 
0  100  0 
0  100  100 
0  100  200 
100  0  0 
100  0  100 
100  0  200 


Sensor  Array  2:  0  0  0 

0  0  100 
0  100  0 
100  0  0 
0  100  100 
100  100  0 
100  0  100 
100  100  100 
50  50  50 


^  =  1.0  meters,  100  trials  per  \IC 


run 


SI  Method  Source  Location  RMS  Error, 
Cramer-Rao  Lower  Round 


LOO  37  37  0.99 


SI  Method  Source  Direction  Cosines  RMS  Error 
Cramer-Rao  Lower  Bound 


1.81  8.0  6.7  1.19  6.7  6.0  1.11 


Table  4 

SI  Method  Source  Range  RMS  Error, 
Cramer-Rao  Lower  Round 


RMS  Error,  Cramer-Rao  Lower  Bound  (meters) 
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TRACK  PARAMETER  ESTIMATION  FROM  MULTIPATH  DELAY  INFORMATION 


J.  S.  Abel  and  Khosrow  Lashkari 

Systems  Control  Technology,  Inc. 
1801  Page  Mi  11  Road 
Palo  Alto,  California  94304 


Abstract 

It  is  desired  to  track  the  location  of  an  underwater  acoustic  source  with 
range  difference  measurements  from  a  stationary  passive  array.  Many  times, 
the  array  has  only  one  or  two  sensors,  and  the  multipath  and  intersensor  range 
difference  measurements  are  insufficient  to  localize  and  track  a  source  moving 
along  an  arbitrary  path  [1],  Here,  we  propose  to  track  sources  with  1-  or  2- 
sensor  stationary  passive  arrays  by  making  the  simplifying  assumption  that  the 
source's  patn  can  be  described  by  a  small  set  of  so-called  track  parameters 
and  using  tne  range  difference  information  to  estimate  the  track  parameter  set 
rather  than  the  source  location  as  a  function  of  time. 

In  this  paper,  we  choose  the  track  parameters  to  specify  a  straight-line, 
constant-velocity,  constant-depth  path.  Cramer-Rao  bounds  are  presented  for 
estimating  these  track  parameters  from  the  time  history  of  multipath  and 
intsrsensor  range  difference  measurements.  It  is  shown  that  this  track 
parameter  set  cannot  be  accurately  estimated  from  the  time  history  of  a  single 
multipath  range  difference  without  side  information  (an  independent  velocity 
estimate,  for  instance),  although  multipath  and  intersensor  range  difference 
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measurments  from  a  two-sensor  array  are  generally  sufficient  to  estimate  the 
track  paremeter  set.  Computationally  efficient  techniques  are  presented 
which  estimate  track  parameters  from  range  difference  measurements  taken  from 
1-  and  2-sensor  arrays. 

Monte-Carlo  simulations  are  presented  which  show  that  these  techniques 
have  sample  mean  square  error  approximately  equal  to  the  Crammer-Rao  bound 
when  a  single  multipath  range  difference  and  an  independent  velocity  estimate 
are  available.  The  sample  mean  square  error  is  shown  to  be  in  the  range  of 
two  to  ten  times  the  corresponding  Crammer-Rao  bounds  when  these  techniques 
are  applied  to  2-sensor  range  difference  data. 
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1.  INTRODUCTION 


The  tracking  of  a  radiating  source  by  a  sensor  array  is  a  basic  problem 
in  underwater  acoustics.  In  this  paper,  we  describe  methods  for  tracking 
sources  using  range  difference  measurements  from  passive  stationary  sensor 
arrays  when  the  number  of  sensors  is  small  and  multipath  surface  reflections 
are  present. 

The  signals  recieved  by  the  sensors  of  an  array  exhibit  relative  time 
delays  correspond!' ng  to  a  sources'  location.  The  relative  time  delay  between 
two  sensors,  termed  time  difference  of  arrival  (TDOA)  is  proportional  to  the 
source-sensor  range  difference  (RD),  and  can  be  used  to  locate  a  source  on  a 
hyperboloid  of  revolution  about  a  line  drawn  through  both  sensors  (a  constant 
velocity  medium  is  assumed).  In  a  p-dimensi onal  space,  N  >  p  independent 
RO's  are  required  to  localize  a  source  without  additional  information  [1], 
Accordingly,  a  target  moving  along  an  arbitrary  path  can  be  tracked  in  p 
dimensions  with  N  >  p  RD's  by  smoothing  source  location  estimates  obtained  at 
various  times. 

Much  work  has  been  done  on  the  problem  of  estimating  TDOA's  (a  special 
issue  on  the  topic  is  [2]  and  recent  studies  on  the  performance  of  time  delay 
estimators  include  [3-10]).  However,  there  is  relatively  little  work  on  the 
problems  of  estimating  source  locations  from  TDOA  measurements  [11-15]  and 
tracking  moving  sources  using  TDOA  measurements  from  stationary  passive  arrays 
[16-19]. 

RD  information  from  1-  or  2-sensor  stationary  passive  arrays  is 
insufficient  to  localize  and  therefore  track  a  source  moving  along  an 
aroitrary  path.  However,  if  it  is  assumed  that  the  source  travels  along  a 
path  described  by  a  parameter  set,  the  problem  of  estimating  a  source's 
location  at  many  points  in  time  (i.e.,  the  problem  cf  tracking  a  source)  is 
reduced  to  the  problem  of  estimating  the  parameter  se*-  describing  the  source's 
path  and  may  be  possible  with  relatively  little  RD  information. 

Many  underwater  sources  of  interest  travel  along  strai ght-1 i ne,  constant- 
depth,  constant-velocity  paths,  specified  by  so-called  track  parameters  -- 
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velocity,  depth,  bearing  angle,  and  the  radius  and  time  of  closest  approach  to 
the  sensor  array.  In  this  paper,  computationally  efficient  techniques  for 
estimating  a  sources'  track  parameters  from  RO  measurements  taken  from  1-  or 
2-sensor  stationary  passive  arrays  in  the  presence  of  a  surface  multipath 
reflection  are  presented  and  evaluated. 

The  structure  of  this  paper  is  as  follows.  Section  two  describes  the 
track  parameter  estimation  problem  and  reviews  basic  RD  relationships.  The 
Cramer-Rao  lower  bound  on  the  variance  of  estimating  the  track  parameters 
given  RD  and  side  information  is  derived  in  Section  three.  Track  parameter 
estimation  methods  are  developed  in  Section  four.  Section  five  reports 
computer  simulation  results.  Section  six  contains  concluding  remarks. 
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2.  The  Track  Parameter  Estimation  Problem 

We  consider  the  following  scenario  illustrated  in  figure  1.  There  is  a 
stationary  passive  1-  or  2-sensor  array  listening  to  a  target  in  the  presence 
of  a  surface  multipath  reflection.  The  source  moves  by  the  array  at  constant 
velocity  and  constant  depth  along  a  straight  line.  Range  difference 
measurements  (between  the  direct  paths  and/or  multi  paths)  are  available  at 
various  points  during  some  interval  of  time.  The  track  parameter  estimation 
problem  is  to  use  the  range  difference  measurements  to  estimate  the  parameters 
describing  the  target's  path.  3e1ow,  track  parameters  are  defined  and  RO 
relationships  are  described  for  three  arrays:  single  sensor,  two-sensor 
vertical,  and  two-sensor  horizontal. 

Denote  L  .  and  L  .  as  vectors  containing  the  (x,y,z)  coordinates  of  the 

— Si  St 

sensor  locations  for  sensors  1  and  2.  Denote  £  as  the  vector  of  track 
parameters,  and  define 


9  s  !  v  v  X-r  V-r*  Z ' 

A  <■  x  y  T  J  i  T  ; 


Here,  v  and  v  are  the  x-axis  and  y-axis  source  velocities,  z_  is  the 
x  y  i 

source  depth  and  x^.  and  y  determine  the  time  and  range  of  closest  approach 
of  the  source  to  the  sensor  array.  At  any  time  t.  ,  the  source  location  may 
be  given  as 


(2.2) 


lT  =  0  -t.  I-  p  =  fx_-v  t.  yT-v  t.  zT ' 

-T  l  3  L  T  x  1  J  T  y  l  T  ■ 

0  U 


where  I3  is  the  3x3  identity  matrix.  The  range  difference  between  the 
source  aid  the  sensors  at  time  t-j  is  given  by 

—12^  ~  ”  ^t!2 ' 1  ^  (2.3) 

where  _x(  )  denotes  the  i-th  component  of  t me  vector  x_  ,  tne  vector 

d,,  contains  the  range  differences  for  times  t.  ,  i  =  l,...N,  and  d  .  (i)  is 
— l  ~dj 

given  by 
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.Lt(1)  -L^. 


(2.4) 


Single  Sensor 

In  the  case  of  a  single  sensor  in  the  presence  of  a  multipath  reflection, 
the  track  parameter  vector,  and  the  source  and  sensor  locations  can  be  given 
as  (see  figure  2) 

£  -  [v  XT  yT  zt]T 

LT(i)  =  [xrvti  yT  zt]T 

W  ■  [o  0  2s]T  (2.5) 

Here  the  x-y  plane  of  the  coordinate  system  is  the  ocean  surface,  positive 
depth  is  measured  into  the  ocean,  and  the  coordinate  system  origin  is  at  the 
ocean  surface  above  the  sensor.  Note,  due  to  the  radial  symmetry  of  the 
array,  there  are  only  four  elements  of  £  .  The  RD  between  the  direct  path 
and  multipath  for  the  sensor  location  above  is  equivalent  to  the  RO  between 
the  direct  paths  to  the  sensor  locations  (see  figure  2) 
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(true  sensor) 

ts™  ■  [°  o  -«,r 

(virtual  sensor) 

The  RO  vector  is  given  as 


d  =  <L  "  (2.6) 

—  — m  ~ d 

where  the  i-th  elements  of  d^  and  d^  ,  the  direct  and  multipath  ranges,  are 
given  by 
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d  (i )  3  )L  -  l_( i ) a 
th  — sm  l 
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Two-sensor  Vertical  Array 

Here,  the  track  parameter  vector  and  the  source  and  sensnr  irw-*t-?nc  ;»r 
be  given  as  (see  figure  3) 


£  a 

[v  xT  yT  zT] 

lT(i) 

3  [xrvti  yT  2ri 

$ 

,t 

^sl  = 

t°  0  zsi! 

,N 

V* 

> 

T 

-s2  = 

[°  0  zs2 1  ’  zs2 

(2.8) 


Again,  due  to  the  radial  symmetry  of  the  array,  there  are  only  four  nonzero 
elements  of  o_  .  The  intersensor  RO  vector  is  given  as 


-12  =  nil  "  ^d2 


(2.9) 


where  the  i-th  elements  of  d,.  and  d^n  ,  the  direct  path  ranges,  are  given  by 

— ul  *~u£ 

V”  ‘  "tI  -  t-T(i)' 

,,  .2  2  ,  .2,1/2 

-  [(XT  -  Vti)  +  yT  +  (ZT  -  zsl)  ] 


dd2(i)  "  ,lks2  '  W(ih 

=  Ot  -  vti)2  +  yj  *  fzT 


,2  A/2 


(2.10) 


Note  that  this  intersensor  RD  is  the  same  as  the  multipath  RD  measured  from  a 
single  sensor  at  ''O.O.jU^  -  z$^)j  for  the  same  set  of  target  track 
parameters  with  p(  a)  given  by  2_  -  z  0i- z  ,)  rather  than  z_  . 
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Two-sensor  Horizontal  Array 

With  a  horizontal  array  there  is  no  longer  radial  symmetry,  and  the  track 
parameter  vector  and  the  source  and  sensor  locations  can  be  given  as  (see 
figure  4) 


1  [VX  vy  XT  yT 


LT(i)  «  [xT-vxti  yT-vyti 

ksi  3  [xs  0  zs]T 

W  =  K  0  zslT 


(2.11) 


The  intersensor  RD  vector  is  given  as 


(2.12) 


where  the  ith  elements  of  d^  and  »  the  direct  path  ranges,  are  given  by 


v1}  3  ^1  ■  W(i)" 


tCxT-vxti-xs)2  +  +  Czrzs)2l1/2 


dd2^  =  :l^s2  ’  — 1 


(2.13) 


Note  that  the  array  is  symmetric  about  the  x-axis,  also  the  intersensor  RD  is 

insensitive  to  the  sign  of  'zT-z  )  . 
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3.  80UN0S  ON  THE  PERFORMANCE  OF  TRACK  PARAMETER  ESTIMATORS 


In  this  section  Cramer-Rao  lower  bounds  on  the  variance  of  estimating 
track  parameters  from  RO  measurements  are  presented.  It  is  shown  here  that 
the  minimum  variance  in  estimating  the  track  parameters  £  =  [v  Xy  Zy]"'" 
from  the  time  history  of  a  single  multipath  RD  is  typically  large  compared  to 
the  RD  variance  unless  side  information  is  available.  It  is  also  shown  that 
side  information  may  not  be  needed  when  a  two  sensor  array  is  used  and 
intersensor  RD  information  is  available. 

Lower  bounds  on  the  variance  in  estimating  the  track  parameters 
£  =  [v  x,  y-j.  Zy]T  are  presented  for  the  case  of  known  velocity  and  a  single 
multipath  RD  available  at  N  points  in  time.  Lower  bounds  on  the  variance  in 
estimating  the  track  parameters  £  =  [v  Xy  yy  Zy]"1"  are  presented  for  the  case 
of  multipath  and  intersensor  RD's  available  from  a  two-sensor  vertical  array 
at  N  points  in  time.  Finally,  lower  bounds  on  the  variance  in  estimating  the 
track  parameters  £  *  [v^  v^  Xy  yT  Zy]T  are  presented  for  the  case  of  the 
multipath  and  intersensor  RD's  available  from  a  two-sensor  horizontal  array  at 
N  points  in  time. 

The  variance  of  an  unbiased  estimator  of  a  parameter  £  is  bounded  below 
by  the  Cramer-Rao  lower  bound  [20-23]: 


var(£)  =  E(s  /)  -  E(l)E(l]T  >  i’1  (3.1) 

where  I  is  the  Fisher  information  matrix  for  the  parameter  £  given  by 

[20]  1 

L  4  io9f(i!*)]}  <3-2> 

where  f(£|]0  is  the  conditional  probability  density  function  of  e  given  x  , 
and  the  estimate  £  is  based  on  data  x  .  If  the  distribution  on  >  is 
parameteri zed  by  z  (e.g.  £  =  g(£)j  ,  tnen  the  Fisher  information  matrix  for 
t  can  be  written  as  [20] 


Finally,  if  £  and  are  drawn  from  independent  distributions,  then  [20] 


I  =1+1 

x..£  £  1. 


(3.4) 


3.1  Single  Sensor  Case 


From  Section  2  and  (3.1),  the  Cramer-Rac  lower  bound  for  estimates  of 
£  =  [v  Xy  y-j-  Zy]T  based  on  N  single  sensor  RD  measurements  is 


where 


var(£.)  >  I 


3d  T  3d 

\  =  Ca£)T  idClE) 


(3.5) 


(3.6) 


where  £  is  a  vector  composed  of  RD  reservations  at  times  t,.  ,  i=l,...,N. 
The  derivative,  ad/3£  is  easily  evaluated  from  (1),  and  is  given  by 


d  ti(xT'vti)Pi  "(xr"vti )pi  "yTpl  "ZTP1  +  zspl 


tN^XT'VtN'°N  ^XT-VtN^pN  ’yTpN  -ZTPN  +  Zspl 


-  1  1  +  1  1 

pi  3  d  ( i )  '  d  (i )  ’  pi  =  d  ( i )  *  d  ( i ) 

-m  — d  -m  — d 


(3.7) 


where 


(3.8) 


And,  assuming  the  RO  estimates  are  unbiased  and  Gaussian  distributed  with 
correlation  matrix  the  information  in  d_  is  given  by  [20] 


I*  “  Rh1 
a  d 


(3.9) 


3.1.1  Ambiguity  Issues 


As  shown  below,  the  Cramer-Rao  lower  bound  in  (3.5)  is  typically  very 
large  comoared  to  the  RD  estimate  variant;  therefore  low  variance  unbiased 


s  V 


■VTJ  s’*/'.'  IJk!'*.  KV.**1  r.'.'*  Vd  t' s's 


track  parameter  estimates  cannot  be  obtained  without  additional  information. 

For  many  sets  of  track  parameters,  during  much  of  the  observation  time, 
2 

the  term  (xT  -  vt.)  will  dominate  d.  ( i )  and  d  (i)  ,  and  1 / d,  ( i )  and 
i  i  - d  — m  — d 

l/dm( i )  can  be  approximated  by, 


d^TT  |xT-vti 


(xT  -vti  )‘ 


v11 


(xJ-vti ) 


(3.10) 


where  k^  and  k^  are  constants  depending  on  the  track  parameters  and  the.^ 
sensor  location.  Using  (3.10)  in  (3.8),  we  see  that  the  four  columns  of 

are  approximately  dependent  on  the  three  vectors. 


|xT-vti 


I  Xj-Vti 


1  xT“vtl I 


(3.11) 


xT-vtN 


|xT-vtN| 


|xT-vtN| 


Therefore,  the  matrix  I  =  f — V  i  ■' — ■  will  be  nearly  singular  and  the 

£ 

resulting  bounds  on  the  minimum  variance  of  track  parameter  estimators  will  be 
large  compared  to  the  variance  in  estimating  the  RD’s. 

This  result  can  be  seen  by  examining  the  target/sensor  geometry.  As 
illustrated  in  Figure  b,  there  are  many  sets  of  parameters  which  result  in 
roughly  the  same  measurements  £  .  Away  from  the  sensor  ar^ay,  the 
hyperboloids  of  constant  RD  can  be  approximated  as  cones.  As  a  result,  the 
differences  in  £  between  tne  case  of  a  target  moving  slowly  close  by  a 
sensor  and  the  case  of  a  target  moving  more  quickly,  further  from  the  sensor 
can  be  small  compared  to  the  deviation  in  estimating  £  .  Consequently,  there 


-V-V-'v-' 


can  be  many  sets  of  'equivalent'  track  parameters  for  a  given  £  and  the 
variance  in  estimating  the  track  parameter  set  from  measurements  of  £  alone 
is  high. 


3.1.2  a  priori  Estimate  of  v 


Many  times,  an  independent  estimate  of  the  target  velocity  is  available, 

through  doppler  measurements,  for  instance  [24].  With  an  a  priori  estimate  of 

the  target  velocity  (an  unbiased  Gaussian-di stributed  estimate  with  standard 
2 

deviation  ^  Is  assumed)  the  Cramer-Rao  lower  bound  on  the  variance  of 
estimating  the  track  parameters  is  given  by 


var(£ 


'  J. 
'T  ZT- 


I  v  xT  yT  zT  |  )  >  1 


•1 

£ 


(3.12) 


where 


3l  r  3d.  T 

I  =  f - 1  T  ( - 1  4-  1  1  I 

£  U,£J  d>3£'  A  3£'  V  3£' 

=  ,'-^T  R-lA  +  _1 -  fjLvJ,JV, 

;3£j  *£  Ud/  2  '•jjjJ  ^3£; 

?v 


;  3 .13) 


where 


il 

3£ 


a£ 

=  [1  0  0  U]  and  -r-  is  given  above. 


When  g ^  is  small  enough,  the  columns  of  I  are  no  longer  linearly 
dependent  and  estimates  of  the  remaining  track  parameters  will  nave  variance 
comparable  to  the  RD  and  velocity  estimate  variances  (depending  on  the  amount 
of  RO  data  available  and  the  track  parameters,  see  section  5,  simulations). 
As  a  result,  the  single-sensor  track  parameter  estimation  method  described 
herein  assumes  an  independent  estimate  of  the  velocity. 


3.2  Two-sensor  Vertical  Array 

When  both  inter-sensor  RO  and  multipath  RD  measurements  are  available 
from  a  two-sensor  vertical  array,  the  Cramer-Rao  lower  bound  for  estimating 

T 

trie  set  of  track  oarameters  o  =  [■/  xT  /-  z-]'  is  given  by 

1*1  i 
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where,  I  is  given  below  using  (3.3)  and  (3.4) 

I  =1  +1  +1 

£  £i  %  %2 

3£i  f  3£i  Sdjp  -r  3d _ o  3d  o  3£,  p 

-  \  :irJ  *  *  t*1)  V*  1 


(3.16) 


where  d^ ,  £2  and  £12  are  the  two  multipath  and  inter-sensor  RD  vectors 

(assumed  statistical ly  independent),  the  matrices  I.  ,  I.  and  I.  are 

— 1  -2  -12 

I,  =  R'1 
-1  -1 

I,  =  R"1 
^  —2 

I,  =  Kil 

-12  -12 

and  the  derivatives  a£^/3£  and  3d_2/3£  are  given  by  (3.7)  and  Sd.^3]!  is 


ti(xT-vti)p{  -(vvti)pi  -np‘i  ’2Tpi  "  Ary +  i~riT 

«  •  •  “dl 


tN(xT-vtN)pN  -(xT-vtN)pN  -yTpN 


”ZT P1  '  .  J  +  d (N; 


ddl(N)'  “dl2 


(3.17) 


where 


-  .  1  1 
pi  =  d^U)  '  d^2(i) 


Note  here  that,  by  the  arguments  above,  the  derivatives  od^/s^  , 

3d_2/3£  and  3£12/3£  each  have  approximately  linearly  dependent  columns. 

However,  unless  the  difference  zj2  -  is  small  compared  to  the  size  of 

the  array,  the  sum  I  +  I  „  +  I  will  have  linearly  independent  columns, 

7  £l  £2  £12 

and  the  Cramer-Rao  lower  bound  on  the  variance  of  the  track  parameters  will  be 
comparable  to  the  RO  estimate  variance  (depending  on  the  track  parameters  and 
the  amount  of  data  available,  see  section  5,  simulations). 
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3.2  Two-sensor  Horizontal  Array 


When  Doth  inter-sensor  RO  and  multipath  RD  measurements  are  available 
from  a  two-sensor  horizontal  array,  the  Cramer-Rao  lower  bound  for  estimating 
the  set  of  track  parameters  £  =  [vx  v^  Xy  yT  Zy]  T  is  given  by 


var(£j  >  I 


(3.18) 


where,  I  is  given  as 

I  =  I  +  I  +  I 
£  2-i  %  hi 


,  3d.  -r  3d.  3d„  o d o  3d.  *  3d .  ^ 

s  ( — i*'i 1  r  r — +  c — 4L)  r  ■  f — +  * — ti.'j  t  f Lfc-i 

3E.  2*  3£  j  '  3£  y  '  3£  >  -  3£  ‘  d.„'3£  ' 


where  d^ ,  £  ,  and  £  are  the  two  multipath  and  inter-sensor  RD  functions, 


I.  ,  I .  and  I  are 
il  —2  —12 

I,  =  Rd1 
-1  % 

r  =  r"1 
^2  -2 

I,  =  RHL 
-12  -12 

3d,  3d.-  3d ,  ? 

and  the  derivatives  -r— -  ,  and  —  are  given  by 


h^rVi)0!  •(xT-%titxs>i  -(yryi!pi  -iT„j  ♦  zxP* 

•  •  •  •  « 

*  •  •  •  • 

* 

tNtXt"'xtN'Xs!l!S  tN-Vl"'yti  )°N  '^T"'  xVxs  )°N  '  C^T"',ytH  i»N  "*1^  T  ZS»N 


■\L;  :.-v\  KV.  V’."  K<' 


(x-r\Vvpi 

• 

• 

VWih  -iWi^sK 

•  • 

•  • 

•  • 

-^T"Vl)pl 

• 

• 

• 

”ZTP1 

+  Zxpl 

• 

•  • 

t,  (yT*v  t.)p„  "(xt"v  t,,±x  Ip,. 
l^T  y  iJl  ^  T  x  N  s  JMN 

• 

■^yT'VytN^pN 

“ZTPN 

• 

+  zsp 

"  tlxs°l  ti(yT'vytl)pl 

•  ■ 

•  • 

"(XT“Vxtl^pl+XSpl  ■^yT"VytN^pl 
•  *• 

•  • 

-(:T  ‘  zs) 
• 

•  • 

•  • 

^_N^XT"VxtN^pN  "  CNxspN  tN^yT"VytN^pN 

•  • 

'(XT'Vn}pN+XspN  ‘(yT‘ytN^pN 

• 

-(2T-2s)»N 

(3.20) 

Provided  x  is  large  enough,  we  expect  the  sum  I  +  I  ,  +  I  will  not 

-S  ii  n.2  -12 

contain  linearly  dependent  terms  and  the  Cramer-Rao  lower  bound  on  the  track 

parameter  estimate  variance  will  be  comparable  to  the  RD  variance,  (depending 
on  the  track  parameters  and  the  amount  of  data  available,  see  section  5, 
simulations) . 


4.  Track  Parameter  Estimation 


In  this  section  closed-form  equati on-error  methods  are  developed  for 
estimating  track  parmeters  from  a  single  multipath  RD  and  a  priori  knowledge 
of  source  velocity.  These  methods  are  extended  to  estimate  track  parameters 
using  intersensor  RD  measurements  from  2-sensor  vertical  and  horizontal 
arrays.  Finally,  methods  are  presented  for  estimating  track  parameters  from 
intersensor  and  multipath  RD  measurements  from  2-sensor  vertical  and 
horizontal  arrays  without  a  priori  knowledge  the  source  velocity. 

4.1  Single-Sensor  Track  Parameter  Estimation:  The  Equation-Error  Approach 


In  the  case  of  a  single  sensor  in  the  presence  of  a  multipath  reflection, 
the  functional  form  of  the  multipath  RD  is  known  in  terms  of  the  track 

A  A 

parameters.  Therefore,  a  fit  of  a  model  d_  based  on  track  parameters  £  can 
be  made  to  the  measured  d  .  The  parameter  estimates  would  then  be  chosen  as 


£  3  min  J(d,  d) 

£ 


(4.1) 


for  some  cost  function  J. 


The  minimization  (4.1)  is  over  a  cost  function  which  is  typically  non- 
convex  in  £  .  Therefore,  in  general,  computati onaly  expensive  exhaustive 
search  methods  must  be  used  in  finding  a  solution  of  (4.1).  In  this  case  the 
computational  burden  can  be  overcome  by  choosing  the  cost  function  J  so  that 
functions  of  the  track  parameters  £  appear  as  coefficeints  in  a  linear 
least-squares  minimization,  yielding  a  closed-form  solution  for  £  based  on 
the  measured  RD  d_  the  sensor  location  and  an  a  priori  estimate  of  the 

velocity  v. 

4.1.1  The  Equation  Error  Method 


Recal 1 , 


where, 


d  3  d  -  a, 

—  — m  — d 


•••n  m  ks  w 


V1)  =  uxrvti)  +  yT  +  ( zt+zs)  ] 


(4.3a) 


d^(i)  =  [(xT  -  vti)2  h-  y2  +  (Zrzs)2]1/2  (4.3b) 

o 

What  is  desired  is  an  expression  relating  functions  of  d(i)  to  d  (i)  and 

2  ““  “HI 

<L|(i)  ,  so  that  functions  of  the  unknown  track  parameters  appear  as  linear 
coefficients  which  can  be  estimated  using  least-squares  techniques. 
Manipulating  (4.2)  gives  the  following  relationship 

[d2(i)  -  4(i)  -  4(i  )]2  -  44(i)4(i)  =  0  (4*4) 

In  (4.4)  d  (i)  and  d.(i)  are  replaced  with  d  (i)  and  d\(i)  ,  the  values 
— m  — a  — m  — d 

A  A  A  A  A 

predicted  by  the  estimated  track  parameters  o  =  [v  yT  Zj ]  ,  and  as  the 
delays  d.  are  not  known  precisely,  an  equation  error  is  introduced 

[/(i)  -  -  £{i)f  -  4d^(i)d^(i)  =  £i  (4.5) 

We  assume  an  a  priori  estimate  of  the  velocity  (denoted  by  v)  is  available, 

A 

and  assign  v  =  v  .  Using  (4.3)  in  (4.5), 
d.  (i)  -  d.  (i  )(zs  +  t?v  ) 

+  3d2(i  )ti  v{Xy}  -  4d_2 ( i )  { x 2  +  y2  +  z2}  +  16z2{z2}  =  £i  (4.6) 

A  A  A  T 

where  the  remaining  track  parameters  to  be  estimated  are  [Xy  y^  Zy]  , 
d_(i)  is  the  measured  RD  at  time  t.  ,  and  e_.  is  the  corresponding  equation 
error  [25].  Oefining  z  as  the  vector  with  i-th  element  £.  ,  (4.6)  becomes 


dfd)  -  4d2(i)(z2  n2v2) 


d4(N)  -  4d2(N)(z^  +  t2v2) 


(4.7c) 


-Sd^tDtjV  4d2(l) 

•  • 

•  • 

-3d2 (N)t„v  4d2(N) 


-15z2 


(4.7d) 


Defining  the  cost  function  J  as 


i  A  T 

J  -  e  We 


(4.8) 


where  W  is  a  positive-definite  weighting  matrix,  the  set  of  track  parameters 
minimizing  J  is  given  by 


where 


(q2  -  qf  -  q3)1/2 


(4.9) 


£  =  (  ST  W  S  j'1  ST  W  r 

Clearly,  the  track  parameter  estimates  given  in  (4.9)  are  obtained  with  little 
computational  expense  compared  to  global  search  methods. 

The  estimate  p  given  in  (4,9)  will  be  referred  to  as  the  equation-error 
estimate.  The  following  sections  discuss  properties  of  the  equation-error 
estimate  and  present  methods  for  extending  the  equation-error  estimate  to  the 
case  of  intersensor  RO  data  from  2-sensor  vertical  and  horizontal  arrays. 

4.1,2  The  Equation  Error 

Here  we  study  the  cost  function  minimized  ' n  obtaining  the  equat ion-error 


I’WWT 


I  , 


>  -> 

> 

«• 

V 
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estimate.  Using  equations  (4.4)  and  (4.6)  the  i-th  equation  error  can  be 
written  as, 

e(i)  =  [l2(i)  -  d2(i)  ♦  2dm(i)dd(i)]2  -  4o£(i)d2(i)  (4.10) 

A 

where  d.  is  the  measured  RD  value,  d_  is  the  R0  estimated  by  the  track 

*  A 

parameter  estimates,  and  and  d^  are  the  direct  and  multipath  ranges 
given  by  the  track  parameter  estimates.  Rearrainging  (4.10)  gives 

i(i)  -  [2dm(i)dd(i)  +  [d2(i)  -  d2(i))](d2(i)  -  d2(i))  (4.11) 

-2  2 

The  factor  d  (i)  -  d  (i)  is  the  difference  between  the  square  measured  and 

“  A  A 

estimated  RD‘s,  and  should  be  much  smaller  than  d  (i)d^(i)  ,  the  product  of 
the  estimated  direct  and  multipath  ranges.  Therefore,  the  i-th  equation 
error  can  be  apporximated  as 

*0)  '  2dffl(i)dd(i)(d2(i)  -  d2(i)) 

-  w(i)(d(i)  -  d(1 ) )  (4.12) 

where 


w(  i )  =  2dffl(  i  )d_(  i )  (£(  i )  +  d.(  i ) ) 

Note  that  the  cost  function  j  =  eTe  is  made  small  by  selecting  track 
parameters  which  have  iiwj  as  small  as  possible  such  that  d  approximates 
d_  .  The  factor  w(i)  is  made  smaller  by  selecting  track  parameters  which 
place  the  source  closer  to  the  point  (0,0,0),  recall  ^  =  [00  .  Thus, 

it  is  expected  that  minimizing  the  unweighted  equation-error  (minimizing 
J  =  e_T£  )  will  produce  track  parameter  estimates  with  a  bias  in  the 
direction  fxT  yT  zT]  =  [0  0  0]  . 

Minimizing  j  =  J ^  is  equivalent  to  minimizing  a  weighted  i_2  norm  of 
the  difference  of  the  measured  R0  values  and  the  estimated  RD  values. 
Specificly,  the  minimization  J  =  is  equivalent  to  the  minimization 
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where  £  =  d  -  d  and  \  is  a  diagonal  matrix  with  i-th  row/column  entry 

“  2“  ~ 
given  as  w  (i)  . 


★  ★T 

If  it  is  desired  to  minimize  J  =  e  e  ,  the  above  suggests  a 
weighting  matrix  and  an  iterative  procedure  for  doing  so:  solve  (4.7) 
iteratively  with  weighting  matrix 


Wn+1  s  diag{l/w2(i)} 


(4.13) 


where  diag{x_(i)}  denotes  a  diagonal  matrix  with  x_(i)  as  the 
entry,  and  w( i )  is  computed  using  track  parameter  estimates 
solution  of  (4.7)  minimizing  J  =  Wp  _e  .  If  the  process 
error  minimized  is 


i-th  row/ column 
from  the  n-th 
converges,  the 


★  *T  ★ 

J  »  £  £  (4.14) 


★ 

Note  that  minimizing  J  rather  than  J  places  more  relative  weight  on  the  RD 

estimates  taken  when  the  source  is  closer  to  the  sensor.  The  benefit  of 
•* 

minimizing  J  over  J  may  be  evaluated  by  the  studying  the  Fisher  information 
matrix  (3.12)  for  the  relative  importance  of  RD  measurements  as  a  function  of 
source  range. 

4.1.2  Variance  of  the  Equation-Error  Estimate 

In  this  section  we  present  expressions  for  the  variance  of  the  unweignted 
equation-error  track  parameter  estimates  when  RD  and  a  priori  velocity 
information  is  noisy. 


Velocity  Known  Precisely 

We  first  derive  an  expression  for  the  equati on-error  estimate  variance 
when  the  velocity  is  known  precisely  and  the  RD  vector  is  unbiased  with 

variance  R^  .  When  the  RD  standard  deviation  is  small,  the  variance  of  the 

track  parameters  can  be  given  by  [25] 


var  (£ 
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'3d_;  ^d  ^  3d.j 


(4.15) 
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where  3£/s£  car  be  evaluated  by  the  chain  rule: 


3£  3£  3£ 

3£  "( 3£)  ^"3£^ 


(4.16) 


3£ 


From  (4.7),  —  is  given  as 

d£ 


3£ 

3£ 


0 

1 

fl 

Po 


0 

0 


1 


1 


0  0 

where  is  the  i-th  element  of  £  .  Recall  in  the  noiseless  case. 


(4.17) 


S  £  =  r 


(4.18) 


(4.18) 

by  £ 

gives 

3S 
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3£. 

3!  3" 

+  $  — 

3d 

II 

3£ 
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3d 

(sT  s 

)*lsT{- 

3S 

Id  d 

(4.19) 


3  £ 

+  Id  1 


(4.20) 


where 


3r 


js  a  dia9{4d3(i)  +  <$d(i)(v2t2  +  z2)} 
3S 

—  a  *  diag{-15d(i)tivq1  +  8d(i)q2} 

RO's  Known  Precisely 


(4.211 


When  the  R0  s  are  known  precisely  and  the  velocity  estimate  v  is  unoiased 
2 

with  small  variance  0  ,  the  equati on-error  estimate  variance  can  be  derived 

in  a  manner  similar  to  (4.15)-(4.21) : 


r  -r  ^  f 


«-»-■«  3.2.  2  T 

varQi  -  [v  xp  yT  zT])  '  (^Kt— JT 


3£ 

where  -r-  can  be  evaluated  by  the  chain  rule: 

o  V 

3£  3£  3£  2  T 

37  s  Caa)CavD  +  0  0  °J 

3£  3£ 

the  derivative  is  given  in  (4.17)  ana  the  derivative  — 

evaluated  as  follows.  In  the  noiseless  case, 


S  £  =  r 


Differentiating  (4.24)  by  v  gives 


aS  3£ 
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Estimated  RO's  and  Velocity 

Finally,  if  the  estimates  of  £  and  v  are  independent, 
of  the  equation  error  estimate  p  =  (v  y_  z_l  can  be  given 

w  1  I  i  J 

tne  sum  of  the  variances  calculated  above 
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can  be 

(4.24) 


(4.25) 


(4.26) 


(4.27) 


then  the  variance 
approximately  as 


(4.23) 


4.2  Two  Sensor  Track  Parameter  Estimation 


In  this  section  methods  for  adapting  the  equation-error  estimate  to 
inter-sensor  RD  measurements  from  two  sensor  arrays  are  presented.  In 
addition,  methods  for  estimating  track  parameters  from  two-sensor  arrays 
without  a  priori  information  are  developed. 

4.2.1  Track  Parameter  Estimation  From  Inter-sensor  RD  Estimates 


The  modification  of  the  equation-error  estimate  (4.9)  to  make  use  of 
inter-sensor  delay  information  rather  than  multipath  delay  information  is 
straight  forward.  If  the  sensors  are  placed  in  a  vertical  array,  the 
intersensor  RD  is  functionally  equivalent  to  the  single  sensor  multipath  RD, 
and  only  the  target  and  sensor  depths  need  to  be  adjusted: 
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where 
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where  S  and  _r  are  as  defined  in  (4.7)  with  d_  =  d_  the  intersensor  RD 
estimate  and  W  is  a  positive  definite  weighting  matrix. 

In  the  case  of  a  horizontal  array,  an  equation-error  estimate  can  be 
derived  in  a  manner  similiar  to  (4.2)  --  (4.91.  r-om  (2.121, 
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where  d,„  is  the  intersensor  RO ,  d  is  the  target  range  to  sensor  1  and 
— 14  ~ dl 

d  .  is  the  target  range  to  sensor  2. 

— d2 

Substituting  (2.13)  into  equation  (4.30), 

—12 ( 1 '  '  ^-12(i  )Kxr-vxti  -  xs^  r  ^xs ^•xT”vxti  ) 


+  ^T“Vyti)  +  (Zs'Zl)  1 

+  4x  (x  -v  t.)  =  0 
S  T  x  l 1 


(4.31) 


Denoting  the  track  parameters  to  be  estimated  by  £  =  [v^  Xj  yj  Zj] 
d^,,  as  the  measured  inter-sensor  RD,  (4.31)  becomes 
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where  an  equation  error  z_.  has  been  introduced.  The  set  of  track  parameters 
T  1  . 

minimizing  J  =  G  W  e  are  given  by 
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wnere  v  =  fv  +  v,  ,  is  the  a  or^or-1  velocity  estimate  and 
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and  W  is  a  positive  definite  weighting  matrix.  Note  that  due  to  symmetries  in 
the  array,  the  signs  of  y^  and  v^  ,  and  (z^-z^)  cannot  be  determined  from 
RO  information  alone;  however,  the  sign  of  y  v  can  be  estimated,  and  here, 

■*  A  A  A  I  * 

v  >0  and  yT  takes  the  sign  of  y^.v 


4.2  Combining  estimates  from  two  sensors 


As  shown  in  section  3,  when  RO  information  from  more  than  one  sensor  is 
available,  side  information  may  not  be  needed  to  obtian  low  variance  track 
parameter  estimates.  Below  methods  for  estimating  track  parameters  from 
inter-sensor  and  multipath  RO  information  from  two-sensor  vertical  and 
horizontal  arrays  are  given.  These  methods  are  not  optimal  and  do  not  produce 
ML  (minimum  variance,  unbiased)  track  parameter  estimates,  but  provide  a 
computationally  inexpensive  alternative  to  the  nonlinear  optimization  involved 
in  computing  the  ML  track  parameter  estimates.  If  the  resources  are  available 
to  compute  the  ML  estimate  via  an  iterative  nonlinear  minimization,  for 
example,  these  methods  provide  excellent  starting  points. 

With  data  available  from  two  sensors  placed  in  a  vertical  array,  sets  of 
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estimates,  parameterized  by  v  can  be  constructed  using  (4.9)  and  (4.29).  A 
line  search  over  v  can  then  be  performed  for  the  set  of  track  parameters 
minimizing  some  cost  function.  Emperically,  it  was  found  that  minimizing  the 
following  cost  function  gave  relatively  unbiased  low-variance  estimates  for  a 
large  range  of  noise  levels  for  the  track  parameter  sets  used  in  the  Monte- 
Carlo  simulations,  see  Section  5: 

jv  -  (2T12(»)  -  ^T2('))2  <4-34> 

where  z-j.^7)  is  the  depth  estimate  based  on  v  using  the  intersensor  RD 
information,  and  z-p2 ^ v ^  is  tf1e  dePth  estiamte  based  on  v  using  multipath  RD 
information  from  the  deeper  sensor. 
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5.  Simulation  Results 
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This  section  reports  computer  simulation  results  on  the  performance  of 
equation-error  estimates  applied  to  noisy  RD  and  velocity  estimates.  Sample 
bias,  standard  deviation  and  RMS  error,  defined  by 


V, 

V. 


-  I  1 

sample  bias  ^  tt  £  (e.  -  a) 

i=l 

N  n 

sample  standard  deviation  -  [-^  £  9^  -  [  e^2]^2 

sample  RMS  error  -  [-j-j-  £  (9.  -  a)2]1/2  (5.1) 


..  (where  9.  is  the  i-th  sample  estimate  of  the  parameter  9  )were  calculated  by 

v.  1 

yt  averaging  results  of  100-trial  Monte-Carlo  runs  and  were  compared  to 

theoretically  calculated  values  and  Cramer-Rao  lower  bounds.  Simulations 
■V;  were  implemented  in  the  Ctrl-C*  language  on  a  VAX  11/785  computer. 
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Results  are  given  for  sources  moving  along  straight-line  constant- 
velocity  paths  by  1-  and  2-sensor  arrays.  Multipath  and  inter-sensor  RD 
estimates  were  calculated  by  adding  white  Gaussian  noise  to  the  true  RD 
vectors;  velocity  estimates  (used  in  the  1-sensor  simulations)  were  calculated 
by  adding  white  Gaussian  noise  to  the  true  velocity. 

Table  1  describes  the  enviromental  data  (sensor  locations,  track 
parameters  noise  levels,  etc.)  for  each  of  the  Monte-Carlo  runs  (the  source 
trajectories  and  sensor  locations  are  also  shown  in  Figure  6).  Note  that  in 
the  case  of  a  single  sensor  or  vertical  array  positioned  at  tne  origin, 

a.  -  [vx  vy  Xy  yT  Zj]^  =  [5  2  250  900  170 ]T  is  equivalent  to 

£L  =  [v  xT  yT  Zy]T  =  [5.38  556.3  747  .8  170  ]T  I  and  [-3  1  140  400  1 10  ]T 

translates  to  [3.152  3.349  423.7  1 10  ]T  .  Each  run  used  100  points  of 

data:  one  RD  sample  every  10  seconds  from  time  t]_  =  -490  to  t^gg  =  500 

seconds,  and  the  total  observation  time  is  t  =  t^gg  -  t^  =  990  seconds.  So 


’’Ctrl-C  is  a  high-level  matrix  calcul  ator  language  sold  by  SCT ,  Inc. 
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tnat  all  biases  and  deviations  are  measured  in  meters,  the  track  parameter 
velocity  is  replaced  here  by  the  distance  travelled  during  the  observation 
interval  vt. 


Single  sensor 

Table  2  shows  sample  bias,  standard  deviation,  RMS  error,  and  theoretical 
deviation  (4.28)  of  the  unweighted  equation-error  estimate  (4.9)  applied  to 
multipath  RD  estimates  and  velocity  estimates  from  a  single  sensor  (Monte- 
Carlo  runs  1  thru  4).  Monte-Carlo  runs  1  and  2  show  the  case  of  perfectly 
known  RD  estimates  and  noisy  velocity  estimates,  and  runs  3  and  4  show  the 
case  of  noisy  RD  estimates  and  perfectly  known  velocity. 

In  these  cases,  the  equation-error  estimator  is  seen  to  be  essentially 
unbiased  and  to  exhibit  a  standard  deviation  comparable  to  the  RD  and  velocity 
deviations.  Note  that  the  theoretical  deviations  given  by  equation  (4.28) 
agree  very  well  with  the  sample  RMS  errors.  Accordingly,  it  was  noted  that 
the  track  parameter  estimate  standard  deviation  appears  to  increase  linearly 
witn  RD  and  velocity  standard  deviation.  As  expected  (by  the  discussion  of 
the  equation  error,  Section  4.1),  the  sample  bias  of  the  equati on-error 
estimate  is  negative  in  runs  3  and  4,  indicating  that  the  track  parameter 
estimate  pulls  the  estimated  source  location  closer  to  the  origin  than  the 
true  source  location. 

Table  3  shows  sample  RMS  errors  theoretical  deviations  (4.28),  and 
Cramer-Rao  bounds  (3.12)  of  the  equati on-error  estimator  appled  to  RD  and 
velocity  estimates  from  a  single  sensor  (Monte-Carlo  runs  5  thru  12).  The 
sample  bias,  not  shown  here,  was  noted  to  be  small  compared  to  the  sample 
standard  deviation.  The  RMS  error  appears  to  be  slightly  smaller  tnan 
theroeti ca I ly  predicted,  consistent  with  a  possible  small  cross-correl at i on 
between  the  RD  and  velocity  noises  seen  in  these  relatively  small  sample-size 
Monte-Carlo  runs.  The  sample  RMS  errors  (and  theoretical  deviations'1  appear 
to  be  very  close  to  the  Cramer-Rao  bounds  indicating  tnat  the  equati or-er,*br 
estimate  is  using  the  RD  and  velocity  information  in  an  efficient  nanner  *'o r 
these  cases. 
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Two-snesor  Vertical  Array 


Table  4  shows  sample  bias,  sample  standard  deviation,  sample  RMS  error 
and  Cramer-Rao  bounds  (3.15)  for  equation-error  estimates  applied  to  RO 
measurements  from  a  two-sensor  vertical  array;  a  priori  velocity  estimates 
were  not  given.  The  track  parameter  estimates  were  made  from  RD  measurements 
taken  from  the  deepest  sensor  with  velocity  estimates  determined  by  (4.34). 

The  equation-error  estimate  sample  bias  is  seen  to  be  a  strong  (perhaps 
quadratic)  function  of  the  the  RD  noise  deviation,  and  is  comparable  to  the 
sample  standard  deviation  at  the  larger  value  of  additive  RD  noise.  The 
sample  standard  deviation  appears  to  be  a  linear  function  of  the  RO  noise 
deviation,  and  small  compared  to  the  track  parameter  values.  The  equation- 
error  estimate  sample  RMS  error  is  seen  to  be  about  twice  the  Cramer-Rao 
bound. 

Two-sensor  Horizontal  Array 

Taole  5  shows  sample  bias,  sample  standard  deviation,  sample  RMS  error, 
and  Cramer-Rao  bounds  for  the  case  of  the  source  moving  past  a  horizontal 
array.  The  track  parameters  were  estimated  by  averaging  the  translated  track 
parameter  estimates  given  by  (4.9)  and  (4.3)  with  a  velocity  estimated  by 
minimizing  the  cost  function  given  in  (4.35).  Here,  the  tracx  parameter 
estimates  are  essentially  unbiased  and  have  a  sample  standard  deviation  which 
is  small  compared  to  the  track  parameter  value  and  appears  to  increase 
linearly  with  RD  standard  deviation.  The  estimates  have  an  RMS  error  of  about 
five  to  ten  times  the  Cramer-Rao  bound,  depending  on  the  track  parameters. 
Here,  it  appears  that  the  track  parameter  estimation  method  is  not  using  the 
RD  information  efficiently.  However,  depending  on  the  application  this 
estimator  might  nave  acceptable  performance. 


6.  Summary 


In  this  paper,  we  discussed  the  problem  of  tracking  sources  with 
multipath  and  intersensor  range  difference  measurements  form  1-  or  2-sensor 
stationary  passive  arrays.  RD  data  gathered  from  such  arrays  are  insufficient 
to  track  a  source  moving  along  an  arbritrary  path,  and  here,  the  problem  of 
describing  the  source's  location  as  a  function  of  time  was  reduced  to  the 
problem  of  estimating  a  small  set  of  parameters  describing  an  assumed 
straight-line,  constant-velocity,  constant-depth  source  path. 

Cramer-Rao  bounds  were  presented  for  estimating  the  track  parameters  from 
the  time  history  of  multipath  and  intersensor  range  difference  measurements. 

It  is  shown  that  this  track  parameter  set  could  not  be  accurately  estimated 
from  the  time  history  of  a  single  multioath  range  difference  without  side 
information.  However,  multipath  and  intersensor  range  difference  measurements 
from  a  two-sensor  array  were  seen  to  be  sufficient  to  estimate  the  track 
parameter  set  when  the  sensors  are  appropiately  placed  and  enough  RD  data  is 
avai 1 abl e. 

Linear  least-squares  equati on-error  techniques  were  presented  which 
estimate  track  parameters  from  independent  velocity  estimates  and  multipath 
range  difference  measurements  taken  from  a  1-sensor  array.  Analytic 
expressions  were  developed  for  the  variance  of  the  estimate,  and  Monte-Carlo 
simulations  were  presented  wnich  show  that  these  estimators  are  relatively 
unbiased  and  have  sample  RMS  error  which  is  given  accurately  by  the  analytic 
expressions  and  is  approximately  equal  to  the  Cramer-Rao  bound. 

Line-search  methods  were  developed  to  estimate  track  parameters  from  2- 
sensor  RO  data  when  no  independent  velocity  estimates  were  available.  The 
sample  mean  square  error  of  the  track  parameter  estimates  produced  by  these 
methods  was  shown  by  Monte-Carlo  simulations  to  be  in  the  range  of  two  to  ten 
times  the  correspondi ng  Cramer-Rao  bounds.  These  methods  provide  a 
comoutational ly  inexpesive  alternative  to  the  nonlinear  optimization  required 
in  finding  tne  ML  estimate. 
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Figure  captions 


Figure  1  RD  Measurement  Scenario.  A  radiating  source  is  shown  moving  along 
a  constant-depth,  constant-velocity,  straight-line  path.  The 
sensor  is  shown  recieving  signals  from  the  source  both  directly 
and  from  a  surface  multipath  reflection. 


Figure  2  Single  Sensor  RD  Measurement  and  Source  Track.  Figure  2a  shows 

the  sensor  at  L  .  located  on  the  z-axis  and  the  source  at  L_  . 

— sd  — T 

The  direct  path  and  multipath  lengths  are  shown  as  dj  and  d 

dm 

Note  that  the  multipath  length  d^  is  the  same  as  the  direct  path 

length  from  the  source  to  a  'virtual1  sensor  located  at  L 

— sm 

Figure  2b  shows  the  source's  constant-velocity  motion  parallel  the 
x-axis  and  the  sensor  located  at  the  origin  of  the  x-y  plane. 


Figure  3  Vertical  Array  RD  Measurement  and  Source  Track.  Figure  3a  shows 

sensors  at  L  .  and  L  on  the  z-axis  and  the  source  at  L  . 

—si  -s2  — T 

The  direct  path  lengths  from  the  source  to  the  sensors  are  shown 
as  and  d ^  ;  the  multipaths  are  not  shown.  Figure  3b  shows 

the  source's  constant-velocity  motion  parallel  the  x-axis  and  the 
sensors  located  at  the  origin  of  the  x-y  plane. 


Figure  4  Horizontal  Array  RD  Measurement  and  Source  Track.  Figure  4a  shows 

sensors  at  L  .  and  l  _  at  the  same  depth,  and  the  source  at 

.  The  direct  paths  from  the  source  to  the  sensors  are  shown 

and  have  lengths  d  ,  and  d  ;  the  multipaths  are  not  shown, 
dl  d2 

Figure  4b  shows  the  source's  constant-velocity  motion  and  the 
sensors  located  on  the  x-axis. 


Figure  5  Track  Equivalence.  Two  sources  moving  along  different  constant- 
velocity,  constant-depth,  straight-line  paths  are  shown  here  to 
produce  similiar  sets  of  single-sensor  multipath  RD 
measurements.  Source  1  is  shown  moving  along  a  constant-velocity, 
constant-depth,  straignt-1 i ne  path  producing  RD  measurements 
[d^  d^]  at  times  [t ^  t,  t^]  .  Source  2  is  shown  moving  faster 

at  greater  depth  along  a  constant-velocity,  constant-deptn , 


strai ght-1 i ne  path  producing  the  same  RO  measurements  at  the  same 
times. 

Figure  6  Simulation  Source  Traces.  The  x-y  plane  source  tracks  used  in 
the  Monte-Carlo  simulations  are  shown  as  dotted  lines,  each  dot 
representing  the  source  location  at  the  time  of  an  RD 
measurement.  The  sensor  locations  in  the  x-y  plane  are  shown  as 
squares;  the  vertical  array  and  single-sensor  array  are 
represented  by  the  center  square  and  the  horizontal  array  is 
represented  by  the  two  outer  squares. 
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Table  1 

Monte-Carlo  Runs 


Track  parameter  set  #1:  [vx  vyxT  yT  ^=[5  2  250  900  170] 

Track  parameter  set  #2:  [vx  vyxT  y^  Zj.]=[-3  1  140  400  110] 

Single  Sensor  Location:  (0,0,300) 

Vertical  Array  Sensor  Locations:  (0,0,100),  (0,0,300) 

Horizontal  Array  Sensor  Locations:  (200,0,300),  (-200,0,300) 


Table  3b 

Track  Parameter  Theoretical  Deviation, 
Single-Sensor  Array 


Run 

Theoretical  Deviation  (meters) 

# 

<?X 

ay 

Oz 

Gvt 

5 

1.3 

2.0 

0.37 

10 

6 

4.2 

6.2 

0.91 

10 

7 

11 

17 

3.2 

100 

8 

14 

20 

3.8 

100 

9 

0.22 

2.3 

0.38 

10 

10 

1.9 

4.9 

0.69 

10 

11 

0.47 

20 

3.5 

100 

12 

2.2 

23 

3.8 

100 

12 


1.7 


20 


3.5 


100 


Track  Parameter  Estimate 
Bias  and  Standard  Deviation, 


Table  5a 

Track  Parameter  Estimate  Bias, 
Two-Sensor  Horizontal  Array 


Run 

Bias  (meters) 

# 

A 

xT-xT 

yT-yT 

A 

Z-f-Z  T 

(vrvx)t 

(vY-vy)t 

13 

-0.11 

-0.61 

-0.10 

2.9 

• 

o 

14 

0.49 

1.2 

0.35 

13 

7.6 

15 

-0.49 

-1.4 

-0.25 

6.7 

-2.3 

16 

0.26 

1  ? 

0.45 

-14 

5.3 

Table  5b 

Track  Parameter  Standard  Deviation, 
Two-Sensor  Horizontal  Array 


Run 

Standard  Deviation  (meters) 

# 

<7x 

Oy 

<?z 

<JV  t 

X 

13 

1.1 

4.9 

0.82 

24 

9.9 

14 

10 

49 

8.2 

240 

99 

15 

1.3 

3.5 

0.65 

18 

6.2 

16 

14 

40 

7.2 

200 

68 

Table  5c 

Track  Parameter  RMS  Error, 
Two-Sensor  Horizontal  Array 


*  *  , 

_  ! 


Run 

RMS  Error  (meters) 

# 

<?x 

ay 

Gv  t 

X 

Gy  t 

Y  1 

13 

1.1 

4.9 

0.83 

25 

10 

14 

11 

49 

8.2 

240 

99 

15 

1.3 

3.8 

0.70 

19 

6.6 

16 

14 

40 

7.2 

200 

68 

Table  5d 

Track  Parameter  CRB, 
Two-Sensor  Horizontal  Array 


Run 

CRB  Deviation  (meters) 

# 

CTX 

a, 

<?z 

cv 

Gy  t 

Y 

13 

0.24 

0.90 

0.15 

4.6 

1.8 

14 

2.4 

9.0 

1.5 

46 

18 

15 

0.086 

0.38 

0.065 

1.8 

0.53 

16 

0.86 

3.8 

0.65 

18 

5.3 

■».V 


WHITE  PAPER 


A  FLEXIBLE  SONAR  SIGNAL 
PROCESSING  WORKSTATION 


January  1986 


Prepared  by: 

Julius  0.  Smith 
Yair  Barniv 


Approved  by: 

Yair  Barniv 
Manager, 

Adaptive  Systems  Department 


SUMMARY 


A  sonar  signal-processing  workstation  is  described  which  can  be  quickly 
reconfigured  on  site  to  provide  maximum  responsiveness  to  new  and  changing 
threats.  The  areas  of  new  technology  are 

o  A  high-level  programming  environment  tailored  to  digital  signal 
processing  development  and  implementation; 

o  Sophisticated  software  tools  specifically  suited  to  underwater 
detection,  localization,  and  characterization. 

The  workstation  thus  serves  in  two  very  distinct  ways,  first  as  the  ideal 
algorithm  development  environment,  and  secondly  as  implementor  of  "production 
algorithms." 

Algorithms  are  developed  interactively  in  a  very  advanced  signal 
processing  language  by  a  highly  trained  signal  processing  expert.  A 
successful  high-level  description  can  then  be  "compiled"  directly  into 
numerous  low-level  machine-code  modules  which  can  control  a  diverse 
installation  of  signal  processing  hardware  in  real  time. 

The  workstation  is  simultaneously  capable  of  high-speed  development  when 
adjusting  to  new  threats  and  high-speed  performance  when  implementing  the 
signal  processing  tools  thus  developed. 


1.  INTRODUCTION 


Many  systems  have  been  built  to  exploit  the  latest  and  most  effective 
techniques  for  signal  processing  in  the  context  of  underwater  surveillance. 
Software  development  for  these  systems  has  typically  relied  on  the  tools 
provided  by  computer  manufacturers.  Once  a  system  is  developed  for  a 
particular  task,  it  is  normally  not  easy  to  modify  in  the  field  on  a 
fundamental  level,  and  the  user  may  be  forced  to  choose  among  a  variety  of 
menu  selections,  unable  to  contribute  fundamental  enhancements  to  the  system. 
While  appropriate  for  lowly  skilled  operators,  such  systems  cannot  be 
developed  on  site,  thereby  significantly  under-utilizing  the  experience  and 
insight  of  the  highly  trained  operator.  The  system  proposed  in  this  paper 
can  close  the  gulf  between  designer  and  end-user;  fundamentally  new  approaches 
to  target  detection  can  be  installed  literally  overnight. 

This  paper  outlines  a  workstation  concept  one  level  higher  than  the 
typical:  In  addition  to  providing  a  powerful  set  of  tools  for  acoustic 
underwater  surveillance,  we  will  provide  a  system  maximally  effective  for 
developing  future  surveillance  systems.  This  allows  maximum  flexibility  with 
respect  to  countermeasures  which  may  render  existing  systems  ineffective. 
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2.  FEATURES  OF  THE  SIGNAL  PROCESSING  DEVELOPMENT  SYSTEM 


Present-day  software  tools  for  maximally  advanced  development  consist  of 
the  following  features: 

Interactive  high-level  language  (with  compiler) 

Integrated  editor 

Self-documenting  and  self-checking  applications  tools 
2.1  HIGH-LEVEL  SIGNAL  PROCESSING  LANGUAGE 

Modern  software  development  systems  include  an  interactive  programming 
language  which  supports  the  notation  and  needs  of  the  application  to  the 
greatest  possible  extent.  In  the  case  of  underwater  acoustic  signal 
processing,  this  language  would  support  complex  matrix  operations,  polynomial 
manipulations,  controllable  numerical  precision,  and  operations  applied  to  an 
indefinitely  long  sequence,  to  name  a  few.  The  language  must  also  support  the 
development  of  concise  mathematical  notation  as  close  as  possible  to  that 
used  by  a  design  engineer  working  with  paper  and  pencil.  The  nearly 
mathematical  descriptions  manipulated  by  the  design  engineer  are  executed 
interpretively  in  the  development  stage,  allowing  maximally  high-level 
debugging  and  evaluation,  but  the  final  equations  are  automatically  translated 
into  ultra-fast  machine  instructions  by  a  hard  working  compiler. 

Unfortunately,  there  is  no  language  for  signal  processing  which  achieves 
a  satisfactory  percentage  of  the  above  goals.  The  principal  language  in 
current  use  for  signal  processing  is  Fortran,  although  in  the  future,  ADA  will 
be  extensively  used.  While  proposed  changes  to  Fortran  over  the  next  decade 
may  give  some  support  to  treating  arrays  and  even  matrices  as  primitive 
objects  in  the  language,  there  is  no  real  possibility  of  obtaining  full 
support  of  signal  processing  needs  from  either  Fortran  or  ADA  in  the 
foreseeable  future. 


The  best  obtainable  signal  processing  language  to  date  appears  to  be 
CtrlC  (pronounced  "control -see"),  a  software  product  of  Systems  Control 
Technology  (SCT).  CtrlC  is  an  extension  of  MatLab  (for  "matrix  laboratory" ) , 
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which  is  a  numerical  analyst's  desk  calculator  language  invented  by  Cleve 
Moler  at  the  University  of  New  Mexico.  MatLab  expressions  involve  a  single 
data  type:  matrices  of  double  precision  complex  floating-point  numbers. 

MatLab  expressions  are  extremely  concise  for  linear  algebraic  equations.  For 
example,  typing  "x=b\A"  in  MatLab  solves  the  matrix  equation  Ax=b  using 
Householder  reductions  followed  by  back-substitution  (Or.  Moler  is  a  leader  in 
the  field  of  matrix  numerical  analysis  and  was  involved  in  the  development  of 
EISPAK  for  eigen  analysis).  If  Ax=b  is  overdetermined,  then  a  a  numerically 
robust  version  of  the  least  squares  solution  x  =  (A  A)"  A  b  is  produced. 
MatLab  contributed  fundamental  progress  in  the  development  of  a  high-level 
syntax  for  digital  signal  processing  operations. 

CtrlC  is  the  only  extension  of  MatLab  which  supports  user-defined 
functions.  This  extension  is  critical  for  the  use  of  CtrlC  as  a  high  level 
signal  processing  language.  The  design  engineer  can  very  quickly  build  up  a 
large  library  of  functions  written  in  CtrlC.  Compared  with  writing  Fortran, 
the  CtrlC  environment  is  about  one  order  of  magnitude  faster  for  algorithm 
development.  This  advantage  is  due  to  the  high-level  syntax  and  the  fact  that 
it  is  interactive.  To  date,  there  is  not  yet  a  compiler  for  the  CtrlC 
language.  Automatic  translation  from  CtrlC  to  Fortran  is  one  subtask  of  the 
workstation  development  project. 

2.2  INTEGRATED  EDITOR 

In  addition  to  the  language  itself,  a  modern  software  development  system 
provides  an  "integrated  editor"  which  "understands"  the  high-level  language, 
and  contains  tools  fitted  to  the  same  application.  For  example,  a  signal 
processing  language  editor  would  provide  extensive  display  support  for  digital 
signals  and  system  diagrams.  The  editor  should  support  both  text  and  graphics 
as  program  source.  Text  is  used  to  describe  mathematical  functions  in  a 
concise  yet  maximally  clear  manner,  and  block  diagrams  are  used  to  describe 
whole  systems.  (Block  diagram  notation  is  especially  important  when  making 
extensive  use  of  parallel  computation.)  Diagram  blocks  can  include  text  which 
contains,  for  example,  the  mathematical  expression  for  the  function  of  that 
box,  or  the  name  of  a  signal  processing  primitive  and  a  pointer  to  the  list  of 
parameter  values  for  that  instantiation  of  it.  (Signal  processing  primitives 


very  often  have  quite  a  few  parameters  which  must  be  fixed  by  the  designer. 

It  is  normally  necessary  to  effectively  provide  a  table  of  parameters  for  each 
processing  box.  They  are  easily  managed  via  menu-style  interfaces.) 

In  a  well -integrated  editor,  a  high  level  "expression"  (mathematical 
syntax  or  a  piece  of  a  system  diagram)  can  be  "evaluated"  in  the  editor's 
environment  to  produce  high  level  objects  which  can  be  inspected  or  used  in 
subsequent  operations.  During  system  development,  it  is  never  necessary  to 
leave  the  editor.  All  results  of  all  operations  exist  in  forms  which  can  be 
further  manipulated,  displayed,  or  transformed  within  the  editor.  Later,  when 
the  processing  functions  have  been  identified,  and  a  complete  system  has  been 
configured,  the  collected  functions  and/or  diagrams  are  then  compiled  by 
another  editor  command,  and  the  net  result  is  a  high-speed  version  of  that 
capability.  The  compiled  code  is  then  available  as  a  (fast)  primitive  in  the 
editor.  It  is  also  available  for  installation  as  an  operator's  option  in  a 
menu  tree,  or  as  a  permanent  ongoing  real  time  process  (such  as  needed  for 
automatic  detection  of  a  new  threat). 

SCT  is  presently  committed  to  the  development  of  the  block-diagram 
interface  to  CtrlC.  Development  of  an  integrated  editor  environment  "around" 
this  facility  and  the  CtrlC  language  is  a  subtask  of  the  proposed  workstation 
development  project. 

2.3  SIGNAL  PROCESSING  TOOLS 

In  addition  to  the  language  and  an  integrated  editor,  there  should  be  an 
extensive  library  of  sonar  signal  processing  tools  which  support  the  design 
engineer  in  every  conceivable  way.  For  example,  primitive  functions  should 
exist  for  computing  the  FFT,  power  spectral  density,  autocorrelation  function, 
cross-correlation,  Doppler  correlation,  target  localization  from  delay  or 
Doppler  information,  digital  filter  design,  sampling-rate  conversion,  coherent 
and  non-coherent  sinusoidal  line  tracking,  peak  finders,  statistical 
utilities,  signal  generators  for  testing,  image  processing  utilities,  and  so 
on.  These  primitives  introduce  higher  level  data  types  which  can  be  given 
"methods"  describing  how  to  display  them,  or  even  how  to  interpret  them  in 
high-level  expressions. 
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All  of  the  signal  processing  primitives  above  and  more  have  been 
developed  in  Fortran  and  CtrIC  at  SCT.  Due  to  lack  of  a  wide  customer  base, 
only  the  most  general  signal  processing  utilities  are  being  sold  in  the 
commercial  distribution  of  CtrIC.  With  a  relatively  small  level  of  effort,  a 
complete  installation  of  sonar  signal  processing  tools  into  CtrIC  can  be 
carried  out  at  SCT.  This  alone  would  seem  to  provide  the  most  advanced  sonar- 
signal  -processing  detection  workstation  in  the  world. 

Part  of  the  process  of  integrating  signal  processing  tools  into  CtrIC  is 
to  supply  intelligent  defaulting  of  unspecified  arguments,  and  consistency 
checking  for  the  arguments.  Here  the  algorithm  designer  can  save  the  user 
from  simple  mistakes  such  as  undersampling  (which  produces  aliasing)  and  "do 
the  right  thing"  when  various  optional  parameters  are  omitted  in  the 
interactive  calling  sequence.  Hardening  of  the  user  interface  in  this  way  is 
expensive  but  very  valuable.  It  and  the  function  documentation  (described 
below)  will  account  for  nearly  all  of  the  time  it  will  take  to  integrate 
existing  SCT  acoustic  surveillance  tools  into  CtrIC. 

2.4  SELF-DOCUMENTATION 

Finally,  given  a  powerful  language  and  integrated  editor,  together  with  a 
comprehensive  set  of  primitive  signal  processing  functions,  there  must  be  a 
very  extensive  online  documentation  facility.  This  documentation  should  allow 
the  expert  engineer  to  sit  down  and  get  to  work  immediately  without  reading 
any  manuals  and  with  a  minimum  of  "database  queries"  on  his  part.  Fast 
acquaintance  with  the  enormous  library  of  instantly  available  primitives  is 
made  possible  by  "keyword  searching"  through  the  function-names,  or  any  of  the 
various  levels  of  function  documentation  described  below.  If  he  knows  the 
name  of  the  function  he  wants  (or  the  first  few  characters  of  it),  he  can  ask 
for  a  display  of  the  documentation  for  that  function  at  one  or  more  levels. 

He  can  be  in  middle  of  typing  a  call  to  a  function  in  an  expression  and 
request  information  (such  as  a  description  of  the  function's  arguments)  to  be 
displayed  in  a  side  window  without  having  to  leave  the  expression  he  is 
typing. 

Each  primitive  function  should  be  documented  on  several  levels:  (1)  a 


one-line  description  of  what  the  function  does,  (2)  a  list  of  argument 
mnemonics  and  their  default  values,  (3)  a  list  of  one-line  descriptions  of  the 
function  arguments,  (4)  a  full  description  of  the  function's  purpose  and 
useage,  (5)  high-level  source  code  for  the  function,  and  (6)  application  notes 
and  examples.  Cross-references  can  be  generated  to  pre-existing  systems  using 
any  given  function  (to  allow  inspection  of  parameters  for  comparable 
situations  in  the  past). 

2.5  PROJECTED  BENEFIT 

In  conclusion,  development  of  the  high-level  workstation  we  propose 
promises  the  following  benefits: 

o  Subcontractor  signal  processing  development  will  produce  better 
results  per  dollar  by  an  order  of  magnitude  or  more. 

o  The  time  required  to  configure  a  new  signal -processing  system  (made 
necessary,  for  example,  by  the  need  to  track  new  types  of  targets), 
will  be  reduced  by  at  least  an  order  of  magnitude. 

These  benefits  are  achieved  by  providing  the  best  existing  software 
development  tools,  uncompromising  high-speed  implementation,  and  a  complete 
battery  of  all  relevant  signal  processing  techniques. 
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3.  SCT  BACKGROUND 

The  need  for  the  advanced  development  environment  discussed  here  has 
become  increasingly  apparent  after  many  man-years  of  signal  processing 
software  development  at  SCT.  SCT  has  accumulated  an  extensive  capability  in 
signal  processing  and  has  devoted  special  emphasis  over  the  past  several  years 
to  passive  underwater  surveillance  systems.  In  addition  to  a  full  complement 
of  sonar  signal  processing  utilities  (many  implemented  on  an  FPS  array 
processor),  SCT  has  developed  an  extremely  powerful  dynamic  programming 
approach  to  multitarget  tracking.  In  related  projects,  extensive  capabilities 
in  digital  image  processing  have  been  developed  over  the  years. 

The  comprehensive  interactive  programming  environment  will  allow  rapid 
prototyping  of  signal  processing  systems  on  the  basis  of  experiments  with 
realistic  data.  CtrlC,  developed  by  SCT,  provides  an  initial  start  in  this 
direction.  While  users  of  CtrlC  are  extrememly  happy  with  it,  the  next 
implementation  will  be  far  more  powerful. 

In  the  following  two  subsections,  some  existing  signal  processing 
utilities  and  some  CtrlC  capabilities  are  listed. 

3.1  EXISTING  SONAR  SIGNAL  PROCESSING  TOOLS  AT  SCT 

For  the  past  several  years  the  Adaptive  Systems  Department  of  SCT  has  had 
ongoing  contracts  with  the  Navy  and  others  to  develop  new  digital  signal 
processing  techniques,  mainly  in  the  area  of  underwater  surveillance.  The 
emphasis  has  been  on  spectrum  analysis  and  source  localization  from  multiple 
hydrophone  recordings.  As  a  result,  the  following  SOFTWARE  capabilities  now 
exi st: 

-  Power  spectral  density  estimation  using  the  FPS  array-processor  (AP) 

-  Estimation  of  crosscorrelation  and  autocorrel ation  functions  on  the  FPS  AP 

-  SCOT  and  PHAT  normalized  cross- spectra  and  correlate  on  FPS  AP 

-  Synthetic  FFT  on  the  FPS  for  narrowband  "zoom"  (The  FHo  maximum  is  8K) 

-  Arbitrary  FFT  window  generator 

-  General  digital  filtering  on  the  FPS  AP 
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-  Sampling  rate  conversion  on  the  FPS  AP 

-  Digital  filter  design  software  to  support  the  above 

-  Modified  Yule  Walker  "high  resolution"  spectrum  analysis 

-  Doppler  cross-correlation  (a  broadband  tracking  technique) 

-  A  multi  target  localization  program 

-  Adaptive  tracking  of  multipath  delay  and  multipath  reflection  coefficient 

-  Adaptive  delay  estimation 

-  Synthetic  signal  generation  (1  or  2  sinusoids  in  white  or  filtered 
noise),  Gaussian  noise  generation,  chirp  sinusoid,  impulse 

-  Spectral  noise- floor  equalizers 

-  Synthetic  Doppler  generator 

-  Signal  statistics  measurement 

-  Signal  interpolation 

-  The  Constant-Modulus  Algorithm  for  multipath  estimation  and  channel 
equal ization 

-  Recursive  maximum  likelihood  method  for  adaptive  filtering  and  estimation 

-  Overdetermined  Instrumental  Variables  method  for  adaptive  filtering 

-  Maximum  A  Posteriori  Line  Extraction  (MAPLE)  for  coherent  narrowband  case 

-  ADEC  for  tracking  noncoherent  sinusoids  in  noise 

-  IEEE  Programs  for  Digital  Signal  Processing 

-  EISPACK,  LINPACK  matrix  analysis  packages  (e.g.  least  squares, 
eigenvalues) 

-  Graylevel  2D  display  (for  spectrograms,  correlograms,  and  the  like) 

-  Graylevel  output  drivers  for  the  HP  LaserJet,  Printronix,  and  DeAnza. 

-  Calibrated  line- track  display  suitable  for  overlay  on  graylevel  plot 

-  Utilities  for  reading  and  unblocking  sonar  data  received  on  magnetic  tape 

-  Support  utilities  for  data  10  in  several  formats. 

Conversion  programs  exist  for  going  between  any  pair  of  formats  below: 

16-bit  packed  integers  (typical  for  raw  data) 

32-bit  floating-point  (typical  for  intermediate  data) 

32-bit  complex  floating-point  (used  for  some  intermediate  data) 

CtrlC  binary  format  (can  load  any  data  into  CtrlC  for  examination) 
MassLib  format  (SCTLIB  image-processing  tools  require  this  format) 
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All  of  the  above  programs  are  written  in  Fortran  or  CtrlC  and  currently  run  on 
the  Vax  11/785.  In  all,  we  have  several  hundred  more  general-purpose  signal 
processing  utilities  which  combine  to  give  a  complete  and  well-rounded  set  of 
tools.  It  is  high  time  for  a  new  phase  aimed  at  integrating  these  tools  in  a 
manner  that  will  put  them  effectively  in  the  hands  of  the  capable  user. 

3.2  SUMMARY  OF  CURRENT  CtrlC  FEATURES 

CtrlC  is  currently  aimed  primarily  at  the  control  design  and  general 
signal  processing  community.  A  feel  for  the  present  CtrlC  environment  is 
gleaned  by  typing  HELP: 

S  DO  CTRLC 
[>  help 

Ctrl-C  Help  Facility 

Starting  &  Ending  Examples  Commands  Summaries 

Basics  Graphics  &  Plotting  Procedures  &  Functions 

Files  &  Data  Matrix  Computations  Control  Systems 

Filtering  &  Estimation  Modeling  &  Simulation  Signal  Processing 

Identification  Numerics  &  Utilities  Personalizing  Ctrl-C 

o  Type  help  and  the  first  unique  letters  of  the  topic  for  info 
o  Type  help  and  the  first  letters  of  AMY  topic  to  check  for  info 
o  If  you  are  new  to  Ctrl-C,  type  'help  intro' 
o  For  more  information  on  how  to  use  the  Help  Facility,  type 
‘help  help* 

o  More  detailed  information  is  available  in  your  Ctrl-C  manual 

Appendix  A  contains  a  complete  list  of  the  functions  presently  delivered  with 
CtrlC. 
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4.  AN  EXAMPLE  SCENARIO 


In  this  section  we  try  to  convey  a  "feel1'  for  the  proposed  system  by 
describing  an  imagined  scenario  in  which  an  outside  consultant  has  been  called 
in  to  develop  a  new  processing  technique  using  the  proposed  signal  processing 
development  system.  To  bring  out  as  many  novel  aspects  of  the  system  as 
possible,  we  assume  the  consultant  has  no  prior  experience  with  it.  He  has 
been  told  only  that  he  can  push  the  'help'  button  if  he  has  any  questions 
about  the  system  or  what  it  can  do. 

The  Scenario 


An  underwater  acoustics  expert  has  been  asked  to  take  a  look  at  some 
unidentified  underwater  signals  which  have  been  recorded  by  a  sensor  array. 

It  is  thought  that  these  sounds  are  emanating  from  a  new  type  of  underwater 
vehicle.  The  objectives  are  to  infer  the  power  spectrum  of  this  new  source 
(as  a  function  of  its  orientation,  range,  velocity,  and  acceleration),  and  to 
develop  software  for  detection,  localization,  and  tracking  in  real  time. 

The  consultant  sits  down  at  the  console  and  presses  the  help'  key.  An 
"editor  window"  is  formed  containing  a  summary  of  further  help  topics.  (This 
is  characteristic  of  the  way  all  requests  for  information  are  handled:  a  new 
editor  window  is  created  initially  filled  with  the  requested  information.) 
The  help  topic  summary  might  look  as  follows: 

Select  any  of  the  following: 

*  INTRO  -  Orientation  for  the  new  user 

*  HELP  -  How  to  use  the  HELP  facility 

*  LANGUAGE  -  High  level  signal  processing  language  tutorial /summary 

*  TOOLS  -  Heirarchical  synopsis  of  available  signal  processing  tools 

*  SYSTEM  -  Documentation  of  system  utilities  (editor,  windows,  etc.) 

(Select  by  pointing  to  desired  line  with  mouse  and  pressinq  mouse  button, 
or  type  capitalized  keyword  above  followed  by  the  RETURN  key.) 


The  user  then  "elides  the  mouse"  on  the  IMTRO  line  (or  types  INTRO 
<RETURN>).  The  introduction  provides  a  general  overview  of  the  system  and 
gives  enough  initial  information  to  allow  the  user  to  proceed  with  the  other 
initial  help'  choices. 

The  format  of  each  help  topic  is  hierarchical .  The  goal  is  to  make  it 
possible  for  a  new  user  to  obtain  any  level  of  detail  in  a  logical  manner,  or 
allow  an  experienced  user  to  get  a  terse  summary  of  features  with  one  or  two 
mouse  clicks.  The  menus  at  each  level  of  the  tree  are  arranged  so  that  the 
experienced  user  gets  to  the  most  often  needed  summary  windows  as  fast  as 
possible  (e.g.  the  mouse  does  not  have  to  be  moved).  In  all  cases  the  help 
information  appears  in  a  separate  window  away  from  the  "working  window"  in 
which  the  user  is  developing  software. 

The  first  time  user  will  want  to  peruse  the  documentation  for  at  least  an 
hour  or  so,  skimming  through  much  of  the  tutorial  level  information  and 
studying  the  terse  summaries  in  more  detail.  The  signal  processing  language 
will  take  the  most  study  time,  although  it  will  be  very  natural  to  use  (being 
closer  to  true  mathematical  notation  than  typical  programming  languages  and 
having  less  general  data  structures  and  data  types  to  learn  about).  The 
language  tutorial  will  include  examples  and  "worksheets"  where  the  new  user 
can  practice  language  features.  The  TOOLS  help  subtree  is  organized  into 
categories,  much  as  the  initial  help  display  of  CtrlC  is  currently  presented 
(see  above).  All  told,  it  should  not  take  more  than  one  day  to  become  fully 
familiar  with  the  general  capabilities  of  the  system  and  how  to  find  specific 
information  when  needed. 

So,  to  work.  Our  acoustics  expert  wants  first  to  see  power  spectrum 
estimates  from  the  three  nearest  sensors  over  time  intervals  centered  about 
CPA  (closest  point  of  approach).  He  therefore  types 

?SDect<cr> 

where  <cr>  denotes  the  RETURN  key,  and  a  window  appears  containing  a  list  of 
all  tools  with  the  word  fragment  'spect'  in  their  titles.  Similar  commands 
create  a  list  of  tools  with  'spect'  in  their  one- line  summaries,  argument 
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documentations,  or  full  documentations.  In  this  way  the  user  can  control  the 
depth  of  the  keyword  searching.  In  addition  to  keywords,  he  can  specify  so- 
called  "regular  expressions"  for  the  search  key.  Regular  expressions  are  a 
notation  for  complex  pattern  matching  used  in  sophisticated  text  editors. 

The  function  'spectrogram1  looks  promising,  so  he  clicks  the  mouse  on 
that  line  of  the  help  window,  or  he  simply  types 

spectrogram  ? 

where  ?  denotes  <C0NTR0L>?.  A  function  name  followed  by  ?  is  a  request  for 
a  help  window  on  that  function.  At  any  time  a  quick  summary  of  the  various 
1 ? *  options  can  be  created  by  pressing  the  <HELP>  key,  selecting  HELP,  and 
then  selecting  HELP  SUMMARY.  This  window  can  be  left  open  in  the  upper  right- 
hand  corner  of  the  screen  for  constant  future  reference  if  desired.  The 
screen  is  large  enough  for  several  windows  to  be  open  at  once.  Windows  can  be 
shuffled  around  and  resized  at  will.  They  can  even  be  "buried"  and  later 
"revived"  like  papers  on  a  desk. 

The  next  task  is  to  apply  the  spectrogram  function  to  each  of  the  three 
closest  sensors.  Suppose  the  passive  sonar  data  for  10  channels  exist  in 
files  Cn.DAT,  for  n=l ,2, . . . ,10.  The  time  of  CPA  in  each  channel  is  estimated 
by  typing  (notice  that  the  text  to  the  right  of  //  represents  comments) 

EnergyEr.v(x,span[5])  :=  {  //  Define  new  function 

//  Return  energy  envelope  =  unnormalized  moving  average  of  x(i)**2 
xs=x.*x;  //  Array  of  x(i)**2 

FOR  i=l: LENGTH(x) , . .  //  Lowpass  filter  (LENGTH  is  intrinsic) 

y( i )=SUM( xs( i-span+1 : i ) ) ;  //  xs(i)  =  0  for  i<l.  (SUM  is  intrinsic) 

RET'JRN(y);  //  Signal  envelope  estimate 

} 

FOR  i =1 : 10 ,  [e( i ) ,cpa( i ) ]=MAX(EnergyEnv(X( : ,i )=GetSig( "C"&i2s( i ) ) ) ) ; 

The  above  is  equivalent  (except  for  the  optimal  argument  'span'  whose  default 
value  is  5)  to  the  following  Fortran: 


SUBROUTINE  ENERGYENV ( Y ,  X ,  N ,  I  SPAN .WORK ) 

DIMENSION  Y (H )  ,X (N) .WORK (N) 

DO  10  1=1, N 
WORK (I)  =  X ( I  )**2 
CONTINUE 
DO  20  1=1, N 
K  =  I 

SUM=W0RK(K ) 

DO  30  J  =  l, I  SPAN- 1 
K  =  K  -  1 

IF  (K.LT.l)  GO  TO  40 
SUM  =  SUM  +  WORK  (K ) 

CONTINUE 
Y ( I )  =  SUM 
CONTINUE 
RETURN 
END 

PARAMETER  NMAX  =  10000 
DIMENSION  E ( NMAX ),ICPA( NMAX) 

DIMENSION  TEMPI ( NMAX ) ,TEMP2 ( NMAX ) , TEMP3 ( NMAX ) , X ( NMAX , 10 ) 

NCHANS=10 

DO  10  1=1 ,NCHANS 

CALL  GETSIG(N, TEMPI, MAKENAMECC'  , I), NMAX) 

DO  15  J=1,N 
X(J,I)  =  TEMPI  (J) 

CONTINUE 

CALL  ENERGYENV (TEMP2 , TEMPI ,10,5, TEMP3 ) 

EMAX  =  TEMP2{1) 

MAXLOC  =  1 
IF(N.LE.l)  GO  TO  30 
DO  20  J=2,N 

I F ( EMAX . GE . TEMP2 ( J ) )  GO  TO  20 
EMAX  =  TEMP2 (J ) 

MAXLOC  =  J 
CONTINUE 


ICPA(I)  =  MAXLOC 
10  CONTINUE 
END 

Clearly,  Fortran  is  not  nearly  as  expressive  as  the  (slightly  extended) 
CtrlC  syntax  above.  In  addition  to  eliminating  manual  handling  of  loop 
counters  in  many  situations,  the  CtrlC  language  provides  for  any  size  of  input 
whereas  the  Fortran  version  must  predefine  a  maximum  array  size  (and  actually 
allocate  that  much  core  memory  whether  or  not  it  is  ever  used). 

If  a  Fortran  version  is  needed  (e.g.  for  portability  or  compilation  on  a 
fast  processor),  the  CtrlC  code  can  be  automatically  translated  into  Fortran 
by  typing,  for  example, 

Fp  =  FORT(EnergyEnv) ;  //  Fp  is  a  string  containing  the  Fortran 

DISP(Fp) ;  //  Print  the  Fortran  on  the  terminal 

SAVE  Fp  >FortSubs.F77  -append  //  Write  to  disk 

If  the  above  is  being  done  often,  a  simple  function  'SaveF77'  can  be 
quickly  defined  so  that  the  command  'SaveF77(EnergyEnv) '  will  accomplish  the 
same  result. 

The  experienced  user  would  have  the  function  'EnergyEnv'  already  defined, 
and  finding  the  time  of  CPA  would  involve  only  one  line  of  CtrlC.  The  user 
would  have  specified  sometime  earlier  (typically  in  his  automatic  startup 
options)  that  the  file  'FortSubs.Olb'  should  be  searched  for  any  undefined 
functions.  'FortSubs.Olb'  is  library  of  COMPILED  Fortran  modules  with  a 
directory  of  entry  points  at  the  beginning.  The  directory  allows  random 
access  within  the  file  for  fast  selective  loading  as  needed.  After  the 
directory  has  been  searched  the  first  time,  the  directory  is  retained  in  core 
to  speed  up  any  future  searches  for  unrecognized  functions. 


The  three  closest  sensors  are  now  found  by  typing 

key  =  1:10;  //  Passive  sort  key 

[Ed,CPAd,key]  =  SORT(E,CPA,key) ;  //  Ed  =  E  in  descending  order. 

//  Remaining  args  sorted  passively. 

Closest  =  key(l:3);  //  Take  three  closest 

The  channel  numbers  of  the  closest  three  sensors  (as  measured  by  maximum 
signal  energy  at  CPA)  reside  in  the  length  three  array  'Closest1.  A 
spectrogram  of  the  closest  sensor  is  obtained  as  follows: 

SI  =  Spectrogram(X( : ,C1osest(l) ) );  //  Take  numerous  default  parameters 

Show(Sl);  //  Display  spectrogram  (2D  gray  or  3D) 

The  default  spectrogram  parameters  provide  a  standard  "Lofargram"  which 
is  a  smoothed  short-time  power- spectral -densi ty  estimate  (power  density  versus 
time  and  frequency).  The  default  arguments  can  be  overridden  in  either  a 
"sticky"  or  "nonsticky"  fashion.  Sticky  overrides  become  the  defaults  in 
future  calls  to  the  function.  A  list  of  all  of  the  arguments  to  Spectrogram 
and  their  defaults  is  obtained  in  a  side  window  by  typing  "Spectrogram(  ?" 

( <C0NTR0L>?  inside  a  function  call).  The  list  of  defaults  can  be  edited  right 
there  in  the  documentation  window,  and  a  command  exists  which  installs  the 
edited  argument  information  as  new  defaults.  Nonsticky  defaults  are  listed  to 
the  right  of  the  original  defaults.  Sticky  defaults  are  obtained  by  typing 
over  the  old  default. 

Having  looked  at  the  spectrograms,  our  consultant  has  determined  that 
there  is  very  little  narrowband  information  associated  with  the  source,  but 
that  there  is  a  lowpass  broadband  "shelf"  extending  to  about  800  Hz.  By 
comparing  the  measured  spectrograms  against  "noise  cuts"  one  hour  earlier 
(again  quickly  accomplished  interactively  in  a  few  minutes),  he  has  determined 
that  there  is  no  useful  information  above  2KHz.  Therefore,  it  will  be  good  to 
lowpass  filter  and  reduce  the  sampling  rate  of  the  sensor  channels  down  to 
5KHz  or  so: 


for  1*1:10,  XR ( : , i )  =  SrConv(X{ : ,i ) ,5000/Fs) ;  //  Fs  =  original  sampling  rate 
Fs=5000’ 

X=Xr; 

CLEAR(Xr) ;  //  Free  up  memory 

Now  the  CPA  times  will  be  used  to  form  an  estimate  of  the  power 
spectrum.  First,  however,  we  plot  the  spectra  at  CPA  in  an  overlay  to  check 
for  bad  data  (due  to  interference,  low  signal  level,  faulty  sensor,  or  the 
like): 

FOR  i=l: 10,  S(:,i)  =  ABS { FFT ( Hammi ng ( X r ( Cpa ( i)-512:Cpa(i)+511))))"2 

S ( : ,  11 )  =  NoiseFloor;  //  Previously  derived  noise  floor  estimate 

Show(S, 'Overlay' ) ;  //  All  columns  of  matrix  S  plotted  in  overlay  fashion. 

Actually,  it  is  more  efficient  to  say  "S( : ,i )=( t=FFT( . . . ) ) .*C0NJG( t)" 
instead  of  "S(: ,i)=ABS(FFT( ...))  2  ,"  but  that's  a  fine  point.  Note  that  each 
spectrum  can  be  displayed  using  a  different  color  in  the  overlay.  Color 
coding  makes  it  possible  to  follow  individual  spectra  in  the  overlay  more 
easily. 

From  the  spectral  overlay,  it  is  determined  that  only  the  closest  seven 
sensors  have  worthwhile  measurements.  We  wish  to  average  these  to  produce  a 
single  estimate  of  the  source  spectrum.  However,  it  is  not  ideal  to  simply 
average  because  some  sensors  have  better  looking  spectra  than  others.  A 
weight  function  is  formed  from  the  energies  at  CPA  as  follows: 

W  =  Ed( 1:7);  //  Weighting  is  proportional  to  energy  at  CPA 

W  =  W/SUM ( W) ;  //  Normalize 

The  estimate  of  the  power  spectrum  is  given  by  a  weighted  average  of  the 
spectrum  at  CPA  in  each  of  the  seven  sensors  closest  to  the  source: 

Starget  =  S(: ,key(l:7))*W;  //  Yes,  vectors  can  be  used  as  indices 

The  estimated  target  PSD  can  be  hand  edited  in  the  waveform  editor  to 
remove  interference.  Also  it  may  be  possible  to  extrapolate  the  source 
spectrum  in  a  physically  plausible  manner  past  the  frequency  where  it  dips 


below  the  noise  floor.  For  example,  if  the  spectrum  is  rolling  off  at  12dB 
per  octave  when  it  merges  with  the  background  noise,  it  is  reasonable  to 
extrapolate  the  12dB  per  octave  rolloff. 

The  consultant  was  mindful  that  the  above  process  of  estimating  a  power 
spectrum  by  averaging  CPA  spectra  might  occur  again.  Therefore,  near  the 
beginning  of  his  session  he  started  up  a  "history  window,"  and  every  keystroke 
typed  is  being  saved  in  this  window.  The  contents  of  the  history  window  can 
be  written  out  to  a  disk  file  at  any  time  and/or  used  as  a  "batch  file"  later 
on.  At  any  time  he  can  "mouse  over"  to  the  history  window  and  make  edit 
changes  (to  eliminate  a  typo  or  digression,  for  example),  or  he  may  even  grab 
a  set  of  earlier  commands  and  execute  them.  At  all  times  the  full  power  of 
the  editor  is  available  in  every  window. 

Mow  we  have  a  preliminary  estimate  of  the  source  power  spectrum.  Using 
it  and  the  noise  cut,  we  can  estimate  a  Wiener  filter  to  be  used  in  making 
broadband  differential  Doppler  measurements: 

Filt  =  Wiener(Starget,NoiseFloor); 

In  reality  a  separate  Wiener  filter  would  be  desired  for  each  time  step  in  the 
short-time  spectra.  Also  it  may  be  desirable  to  edit  this  weight  function  in 
the  waveform  editor  to  downweight  frequency  bands  of  high  "variance"  in  the 
PSD  estimate.  That  is,  the  Wiener  filter  only  penalizes  closeness  to  the 
noise  floor  while  in  reality  there  should  also  be  a  penalty  for 
nonstationarity  (or  "spuriousness")  over  time.  The  spectral  "reliability 
measure"  can  be  computed  quantitatively  or  drawn  according  to  the  judgment  of 
the  signal  processing  engineer.  Whatever  the  ultimate  weighting,  let  Filt’ 
contain  the  final  spectral  window. 

Next  we  attempt  to  localize  the  target  by  estimating  differential  Doppler 
tracks  in  the  seven  good  sensor  recordings.  Differential  Doppler  is  measured 
by  cross-correlating  the  Weiner  and  confidence-weighted  spectra,  after 
resampling  the  frequency  axis  on  a  log  scale.  The  details  are  similar  to  the 
above  examples  and  will  appear  in  a  forthcoming  final  report  on  a  project  of 
this  nature.  (The  work  was  fully  carried  out  within  CtrlC).  Suffice  it  to 


note  that  development  of  broadband  Doppler  correlation  software  was  greatly 
facilitated  by  the  high-level  interactive  programming  environment  provided  by 
the  present  version  of  CtrlC.  The  planned  extensions  to  CtrlC  will  take  it 
far  beyond  the  present  capabilities. 

After  obtaining  differential  Dopplers  plus,  perhaps,  differential  time- 
of-arrivals  (from  temporal  cross-correlations)  and  multipath  measurements, 
along  with  their  associated  weight  functions,  the  data  are  passed  to  a 
comprehensive  multi  target  tracking  facility  (based  on  dynamic  programming  as 
discussed  in  the  next  section).  The  results  of  localization  can  then  be  used 
to  deconvolve  the  propagation  channel  and  obtain  calibrated  spectral 
estimates  as  a  function  of  orientation  and  all  operating  parameters  of  the 
target. 

The  proposed  system  has  been  described  from  the  point  of  view  of  a  new 
user  trying  to  develop  algorithms  almost  from  scratch.  This  aspect  of  the 
system  has  been  stressed  because  it  is  the  most  unusual  feature:  algorithm 
development  is  supported  at  the  highest  possible  level. 

To  translate  a  finished  algorithm  design  into  a  deployable  form,  it  is 
only  necessary  to  "clean  up"  the  code  in  the  history  window,  translate  various 
CtrlC  functions  into  Fortran  (or  other  suitable  "intermediate  code"),  compile 
the  Fortran  for  the  target  machine  (such  as  an  array  processoror  "super 
computer"  configuration),  and  build  a  block  diagram  of  the  final  processing 
system.  Each  block  in  the  diagram  contains  the  name  of  the  function  that  goes 
into  it.  Block  diagrams  are  not  essential,  but  they  are  convenient 
descriptions  of  complex,  parrallel,  real-time  processing  tasks.  They  also 
serve  better  as  "prints"  for  the  final  system.  Using  "software  probes" 
(controlled  by  the  mouse)  windows  can  be  formed  which  display  the  waveform  at 
any  point  in  the  diagram  (while  it  is  executing  in  real  time),  or  values  can 
be  printed  into  a  side  window  and  further  processed  in  an  unrestricted  manner 
much  as  described  above.  Several  windows  can  be  set  up  to  display  waveforms 
at  several  different  points  within  the  network.  A  waveform  display  window  can 
also  be  allocated  to  automatically  display  the  results  of  each  line  of 
interactive  computation  in  the  CtrlC  language.  At  any  time  the  user  may  move 
over  to  the  waveform  and  interactively  manipulate  it. 
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This  discussion  has  only  scratched  the  surface  of  the  capabilities  of  the 
proposed  system.  Needless  to  say,  a  considerable  effort  will  be  required  to 
fully  implement  the  concept.  However,  the  commercial  availability  of  good 
interactive  programming  environments,  together  with  SCT's  experience  and 
progress  to  date  in  this  area,  render  the  project  feasible. 
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5.  DYNAMIC  PROGRAMMING  APPROACH  TO  MULTITARGET  DETECTION 


In  this  section  we  go  into  more  detail  regarding  our  dynamic  programming 
approach  to  target  detection  and  localization.  Essentially,  sensor 
measurements  are  compared  to  expected  measurements  given  the  target  state,  sea 
state,  propagation  models,  and  clutter  conditions.  While  such  an  approach  is 
clearly  as  powerful  as  any  conceivable  for  target  detection  and  tracking,  it 
demands  extensive  signal  processing  and  modeling  support  such  as  can  be 
provided  by  the  proposed  system. 

5.1  MULTITARGET  DETECTION  IN  A  VERY-LOW-SNR  ENVIRONMENT 

In  our  recent  efforts  on  multi  target  localization  we  have  dealt  with 
cases  where  the  frequency  lines  appearing  in  the  Spectrograms  (or  other 
"Grams")  are  relatively  easy  to  recognize — at  least  by  the  human  eye.  In  such 
moderate- SNR  cases  it  is  possible  to  automate  the  line  detection  process  using 
some  off-the-shelf  1 i ne-tracking  algorithms  such  as  ADEC  or  MAPLE  which  we 
have  done  successfully. 

When  the  SNR  of  the  Spectrogram  data  falls  under  some  critical  minimum, 
the  above  line  tracking  algorithms  collapse.  In  such  cases  we  suggest  to 
apply  a  much  more  powerful  approach  which  has  been  developed  at  SCT  for  the 
problem  of  detecting/ tracking  very  dim  airborne  targets  as  perceived  in  a  set 
of  IR  mosaic  sensor  imagery.  This  proprietary  algorithm  is  based  on  what  is 
referred  to  as  Dynamic  Programming  Algorithm  (DPA)  or  Viterbi  Algorithm,  and 
can  be  shown  to  solve  the  detection/tracking  problem  optimally.  In  the 
following  we  will  describe  the  theory  and  possible  implementation  of  this 
al gori thn. 

Although  the  algorithm  is  very  general,  we  prefer  to  describe  it  in  the 
context  of  our  particular  problem  to  be  defined  as  follows.  Having  a  set  of 
sensors  and  the  various  "Grams"  that  can  be  obtained  from  each  (such  as 
LOFARgram,  Spectrogram,  Correlogram  and  Doppler  Correlogram) ,  detect, 
localize,  and  track  targets  in  the  optimal  way  so  as  to  incorporate  all  the 
available  data  and  side  information  into  the  solution. 


For  simplicity,  let  us  assume  that  a  Spectrogram  output  for  each  sensor 
is  the  only  type  of  measurement  available.  Further,  we  are  looking  for 
submarines  which  are  known  to  emit  a  single  narrowband  frequency  and  to  travel 
at  a  constant  velocity  within  some  given  speed  window.  The  unknowns  are  (a) 
the  number  of  targets  (which  may  also  be  zero)  in  the  ocean  volume  under 
surveillance,  (b)  the  speed,  location,  and  direction  of  each  target. 

Let  us  now  imagine  a  single  target  travelling  from  point 
(x^,  yj»  z^)  To  (x,,,  y^,  z2)  during  the  Spectrogram  time  period,  say  10 
minutes.  For  this  target  we  can  predict  the  Spectrogram  lines  that  should 
appear  in  each  one  of  the  sensors  given  their  (x,y,z)  coordinates.  These 


frequency 


SENSOR  n= 1  SENSOR  n=2  SENSOR  n=3 

Fig.  5.1.  Spectrograms  for  3  Sensors 

lines  may  look  like  the  dashed  lines  shown  in  Figure  5.1.  We  can  thus  always 
calculate  the  expected  noise-free  Spectrograms  for  any  number  and  location  of 
sensors  for  any  given  target  trajectory  and  period  of  time. 

If  we  now  discretize  the  time/frequency  Spectrogram  images  by,  say,  10 
seconds/1  Hz  pixels,  we  will  be  able  to  store  the  above  information  in  a 
triply-indexed  (3  dimensional)  table.  The  table  indices  will  be  n  for  the 
sensor  number  and  (i,j)  for  the  time/frequency  pixel  location;  its  entries 
will  be  the  calculated  power  at  each  pixel.  In  other  words,  we  have  mapped  a 
given  general  target  trajectory  into  a  triply-indexed  table  of  expected 
signals.  The  next  step  would  be  to  also  discretize  the  relevant  ocean  volume 
into  cells  and  map  all  possible  (straight)  target  trajectories 
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x  -  (x^  y^  z^ ,  x2,  y2,  z2)  ,  into  tables  using  the  procedure  described 
above. 

We  can  now  "matched  filter"  (MF)  the  actual  (discretized)  Spectrogram 
data  with  the  prepared  bank  of  filters  in  parallel,  and  use  their  outputs  to 
determine  the  likelihood  of  each  trajectory  (see  Fig.  5.2).  It  is  clear  that, 
if  no  target  is  present,  the  outputs  of  all  MFs  will  consist  of  noise  only, 
and  when  a  target  is  present  the  MF  that  "fits"  the  trajectory  will  show  the 
maximum  output.  This  way  we  will  be  able  to  detect  and  track  simultaneously. 


Actual  measured 
spectrograms 


Figure  5.2.  Filtering  the  spectrograms  data  through  all  Filters 
Each  Filter  "fits"  a  Particular  Straight 
Target  Trajectory  in  the  Ocean. 


The  above-described  conceptual  algorithm  amounts  to  an  intensive 
exhaustive  search  over  all  possible  straight  trajectories.  In  practice,  we 
are  of  course  also  interested  in  non- straight  trajectories.  It  is  easy  to 
realize  that  searching,  even  just  for  the  straight  trajectories,  is  a 
monumental  and  infeasible  computational  task.  There  is,  however,  a  feasible 
way  to  do  the  equivalent  of  this  and  other  exhaustive  search  problems;  the 
"magic"  is  to  use  the  Dynamic  Programming  Algorithm  (DPA)  to  do  the  above 
search  iteratively  stage-by-stage  which  requires  many  orders  of  magnitude  less 
computations  and  memory.  The  way  we  can  apply  the  DPA  here  is  the  following. 

We  define  all  possible  short  (say,  1-minute)  and  straight  trajectory 
segments  as  states,  x  -  (x^,y^,z^,  x^,  y£,  z2^  *  w,iere  the  end  points  of  each 
state  are  given  by  their  discrete  ocean-volume  indices.  For  each  such  state 
we  will  calculate  the  table  of  expected  Spectrogram  discrete  measurements  for 
all  sensors.  This  table  will  of  course  be  much  smaller  than  the  one  required 
to  describe  a  complete  trajectory.  Next,  we  define  the  DPA  stages  to  be  1- 
minute  intervals,  for  example.  A  set  of  k  consecutive  states, 


would  thus  describe  a  complete  trajectory  composed  of  short  straight  segments, 
and,  as  such,  would  also  allow  curved  trajectories  (notice  that  the  index  k  of 
x^  is  now  used  as  a  stage  number  (related  to  time)  and  not  as  a  discrete 
volume  index). 

We  now  denote  the  set  of  single-stage  discrete  Spectrogram  measurements 
from  all  sensors  by  z^  (the  solid  lines  in  Fig.  1)  and  the  sequence  of  such 
sets  for  all  stages  (the  complete  Spectrograms)  by 

•  ^  f  > 

K  "  i  Z1  ’  Z2  ’  '  ”  ’  * 

The  other  required  definitions,  assumptions,  and  notations  are: 

pUjZ^)  =  a  posteriori  probability  density  function  (PDF)  associated 
with  a  candidate  trajectory,  given  all  the  measurements 


w  =  random  forcing  function  at  the  k-th  stage;  w  is  assumed 

l\ 

white  and  independent  of  xQ. 

=  measurement  noise,  v  is  white  and  independent  of  w  and 

V 

Using  the  above  notation,  we  can  formulate  the  general  optimum  estimation 
problem  as  follows  [1][2].  Given  a  known  difference  equation  describing  the 
dynamical  system 

xk+l  =  f^xk’  'V  ’ 

the  a  priori  PDF  of  the  initial  states  p'(xQ)  »  the  statistics  of  v,  p(v^ )  , 
and  w,  p(w. ),  i=l,...,k  ,  the  measurements  z^ ,  i=l,...,k  ,  and  the  known 
measurements  relationship 

\  =  h(xk,  vk)  ,  (5.2) 

find  the  sequence  of  states  (trajectory) ,  x  ,  that  maximizes  the  conditional 
POF 


p(Xk|Zk)  .  (5.3) 

Notice  that,  in  our  case,  (5.1)  describes  a  simple  straight-line  extension  of 
states  (or  trajectory  segements),  where  wk  is  the  pixels  quantization  noise, 
and  (5.2)  describes  the  multi-sensor  Spectrogram  measurements,  where  v  is 

IN 

the  measurements  noise. 


The  basic  idea  of  the  proposed  algorithm  is  to  compute  the  likelihood 
function  for  all  possible  trajectori es,  and  choose  those  having  sufficiently 
high  likelihood  as  candidate  detections.  To  see  more  clearly  the  factors  that 
enter  this  function  let  us  consider  what  happens  when  we  advance  by  one 
stage.  Using  Bayes1  theorem  it  can  be  shown  that 


(see  [l]-[5]}.  The  function  P^+JZ^)  is  used  for  normalization  purposes 
only  and  need  not  be  carried  along. 


p(zk+i!xk+i)  is  the  PDF  of  the  measurements  at  the  k+l-st  stage,  given 
that  a  particular  state,  x  ,  contains  a  target.  This  POF  takes  into 
account  the  intensity  of  the  return  in  the  time/frequency  pixels  associated 
with  the  given  state  and  information  about  target  signature.  The  probability 
will  be  higher  if  the  signature  of  the  return  associated  with  state 
xk+l  looks  a  tar9et»  and  Tower  if  it  does  not  have  typical  target 
characteristics.  This  term  will  be  in  our  case  the  output  of  the 
xk+1  matched  filter. 


p ( xk+i I xk ^  1s  the  POfr  a  tar*9et  being  in  state  xk+1  ,  given  that  it 
was  previously  in  state  xk  .  It  is  here  that  we  incorporate  into  the 
algorithm  information  about  target  dynamics.  Since  a  target  can  change  its 
course,  speed  and  location  only  gradually,  high  probability  will  be  given  to  a 


state  xk+1  consistent  with  previous  target  parameters  (course/speed/ 
location),  and  low  probability  will  be  given  to  states  corresponding  to  large 


changes  in  target  parameters. 


We  now  split  the  required  maximization  of  (5.3)  as  follows: 


max  P(\+llZk+1)  =  nax[max{P(Xk+liZk+i)/  •  (5.5) 

Xk+1  Vl  Xk 

Using  (5.4)  and  defining  the  Merit  Function,  I(xk+1)  ,  as  the  inner  maximum 
on  the  righthand  side  of  (5.5),  we  can  arrive  at  the  recursion  (see  [2]  and 
[13]-[16]) 


I(xk+1^  =  p(zk+lK+l^  maxiP!xk+lixk)1(xk^  (5'6) 

xk 

which  defines  a  first-order  Markov  relationship.  Initially,  I(xQ)  =  p(xQ)  Ts 
the  (equal)  a  priori  PDF  for  all  the  states  at  stage  zero.  From  (5.6)  we  see 
that  for  each  state,  xk+1  ,  there  is  a  certain  xk  that  yields  its  maximal 
I ( xk+^ )  •  This  functional  relationship  can  be  expressed  as 

xk  =  xk(xk+l^  •  (5’7) 
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The  algorithmic  interpretation  of  equations  (5.6)  and  (5.7)  is  the 
lowing. 

a)  Define  a  close  set  of  discrete  states  (1-minute  straight  trajectory 
segments)  or  alphabet  which  are  sufficient  to  fully  relate  the  state 
of  the  "system"  at  any  given  time  (or  stage). 

b)  Define  the  discrete  stages  (time  intervals)  so  that  the  system 
advances  from  one  stage  to  the  next  stage. 

c)  Define  the  discrete  Merit  (or  cost)  Function  for  all  states  as  in 

(5.6) . 

d)  Start  with  any  constant  Merit  Function  at  stage  zero  and  calculate 

(5.6)  for  all  allowed  states  at  the  next  stage. 

e)  At  each  stage  replace  the  old  Merit  Function,  I ( x.  )  ,  by  the  updated 

A  ^ 

one,  Kxk+1)  ,  and  store  the  optimal  x^  of  (5.7). 

f)  The  end  resit  at  the  last  stage  is  the  (single)  overall  Merit 
Function  and  the  sets  of  x^  for  all  stages.  This  is  the  point  at 
which  we  finally  perform  the  (last)  outside  maximization  of  (5.5)  and 
find  the  single  state  having  the  largest  Merit  Function.  This  is  the 
last  state  of  the  optimal  trajectory  in  the  case  where  a  single 
target  is  known  to  exist. 

g)  Mow  go  backward  from  this  optimal  last  state  and  use  the  stored 

A 

x^  data  to  retrieve  the  best  previous  state  leading  into  the  optimal 
xq  .  This  process  is  called  "retrieval"  because,  starting  from  the 
best  last  state,  it  is  possible  to  read  off  the  complete  optimal 
trajectory  back  to  the  first  stage.  In  fact  it  is  possible  to  read 
off  the  whole  trajectory  leading  to  any  state  at  the  last  (or  other) 
stage  even  if  the  state  is  not  the  one  yielding  the  maximum  Merit 
Function.  Thus,  in  our  algorithm  we  do  not  perform  the  above  last 
optimization  and,  instead,  we  retrieve  all  tra^  :tories  whose  final 
Merit  Function  crossed  some  threshold.  This  modification  allows  for 


cases  where  it  is  not  known  in  advance  how  many  actual  trajectories 
(if  at  all)  are  present  in  the  ocean  volume  of  interest. 

The  process  described  above  can  be  further  explained  by  reference  to 
Figures  5.3  to  5.5.  Figure  5.3  shows  how  the  initially  constant  Merit 
Function  evolves  throughout  the  stages  until  it  manifests  a  distinct  peak  (for 
a  single-target  case).  The  last  optimal  state  of  the  trajectory  in  this 

A 

example  is  x,  =  st  ,  ,  where  n  is  the  total  number  of  states  in  the 
k  n— o 

alphabet.  The  "signal"  to  si  delobes  "noise"  can  be  used  to  measure  our 
confidence  in  the  detection  process.  Figure  5.4  shows  how  the  DPA  proceeds 
forward  while  looking  backwards  from  each  state  to  determine  its  optimal 
previous  state.  Figure  5.5  shows  the  retrieval  process  for  a  simple  10-state 
case.  States  #2  and  #9  are  retrieved  here. 

5.2  GENERALIZING  THE  DPA 

The  previous  section  explained  the  DPA  application  to  the  single¬ 
frequency  multi -target  case  through  the  use  of  the  Spectrograms  obtained  from 
several  sensors.  The  main  power  of  the  DPA  is  that  the  very  same  general 
formulation  can  take  into  account  all  the  given  a  priori  data  and 
measurements.  For  example: 

the  term  p(x  )  can  relate  any  a  oriori  statistical  knowledge  about 
o 

the  initial  targets'  parameters, 

the  choice  of  the  alphabet  of  states,  x,  can  incorporate  any  known 
location  and  direction  bounds  in  addition  to  the  already  accounted- 
for  speed  bounds, 

the  pUjJx^)  term  (the  output  of  the  multi-sensor  Spectrograms 
matched-filters)  can  incorporate  other  known  frequencies  or  broadband 
sources  as  well  as  Correlogram  and  Doppler  Spectrograms  measurements. 

In  summary,  the  optimality  of  this  algorithm  (assuming  Normal 
measurements  noise)  results  from  the  application  of  matched-filtering  to  the 


data,  and  from  the  incorporation  of  all  the  a  priori  intelligence  or  other 
data  through  the  DPA.  The  DPA  also  makes  this  optimal  detection/tracking 
feasible  in  terms  of  computations /memory  requirements. 

5.3  COMPARISON  WITH  OTHER  TRACKING  ROUTINES 

The  "Maple"  line- tracking  routine  [17]  is  an  example  of  a  "diminished- 
scale"  DPA  application  to  the  line  tracking  problem.  It  uses  matched- 
filtering  like  we  do,  but  it  tracks  frequency  trajectories  defined  in  the 
Spectrogram  image  and  not  physical  submarine  (x,y,z)  trajectories  defined  in 
the  ocean  volume.  Although  this  algorithm  is  optimal  for  what  it  is  doing,  it 
is  far  suboptimal  from  our  viewpoint.  First,  it  does  not  start  frequency-line 
trajectories  everywhere,  and  second,  more  important,  it  does  not  incorporate 
the  signals  from  all  sensors  and  all  "Grams."  Further,  it  cannot  incorporate 
a  priori  knowledge  about  the  expected  trajectory  parameters  or  bounds. 

The  Maple-algorithm  excellent  performance  can,  however,  serve  as  a 
precursor  to  the  performance  which  is  achievable  by  using  our  suggested  full- 
scale  DPA. 
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7.  PROPOSED  STATEMENT  OF  WORK 


7.1  INTERACTIVE  SIGNAL  PROCESSING  WORKSTATION 

The  workstation  project  will  consist  of  the  following  tasks: 
Task  1:  Library  Development 


The  large  collection  of  existing  Fortran  and  CtrlC  signal  processing 
utilities  needs  to  be  organized  into  one  package.  A  uniform  style  of 
documentation  is  needed. 

Task  2:  Module  Interface  Development 


Many  of  the  signal  processing  modules  need  a  "generic  interface"  which  is 
called  interactively  by  the  user.  An  interface  routine  provides  (for  example) 
argument  defaulting,  argument  checking,  dynamic  allocation,  and  automatic  type 
conversions  when  necessary.  Currently,  such  modules  are  used  to  install 
Fortran  subroutines  in  CtrlC. 

Task  3:  CtrlC  Language  Extensions 

Numerous  new  language  facilities  are  planned  for  CtrlC  to  further  enhance 
its  expressive  power.  For  example,  the  macro  capability  will  be  made  much 
more  powerful.  In  addition,  new  high-level  debugging  aids  are  planned. 

Task  4:  CtrlC  Environment  Extensions 


This  task  includes  development  of  all  general  display  and  user  interface 
capabilities.  To  harness  similar  existing  capabilities  as  much  as  possible, 
CtrlC  will  be  ported  to  a  high  quality  interactive  programming  environment 
which  supports  multiple  process  windows.  The  neces$"*y  CtrlC  specific 
enhancements  will  then  be  installed. 


Task  5:  CtrIC  HELP  Extensions 

This  task  involves  installing  the  user-direc  ed  documentation  facility 
described  in  this  paper.  In  essence  the  proposed  help  facility  is  a 
particularly  easy  to  use  relational  database. 

Task  6:  CtrIC  to  Fortran  Compiler 

It  is  desired  to  be  able  to  translate  the  CtrIC  high-level  language  into 
highly  portable  Fortran. 

Task  7:  Block-Diagram  User  Interface 

The  generic  block  diagram  capability  is  being  developed  in  a  separate 
project.  This  task  involves  fully  integrating  the  block  diagramming  facility 
into  the  system  environment. 

Task  8:  The  Integrated  Editor 


Hopefully  the  desired  editor  will  be  largely  provided  with  the  target 
interactive  programming  environment.  Even  so,  some  customization  of  the 
editor  will  be  desirable.  Integration  with  the  block  diagram  facility  is 
desi red. 

7.2  DYNAMIC  PROGRAMMING  FOR  MULTITARGET  DETECTION  AND  TRACKING 

The  DP  algorithm  has  been  developed  at  SCT  for  the  problem  of 
detecting/tracking  aircraft  from  space  as  they  are  observed  in  an  IR  imagery 
set.  The  algorithm  development  got  to  the  point  of  successful  demonstration 
by  simulation,  yielding  results  which  are  in  excellent  agreement  with  the 
theoretical  performance  estimates.  In  the  following  we  suggest  to  develop  a 
similar  algorithm  for  the  ocean  underwater  surveillance  problem  along  the  same 
lines  that  led  us  to  the  existing  algorithm. 


Task  1:  Algorithm  Development 


We  will  develop  the  Dynamic  Programming  Algorithm  for  the  underwater 
surveillance  problem  as  outlined  in  the  preceding  text.  This  will  include 
defining  states,  stages,  transitions  and  measurements  functions,  and  the  Merit 
Function. 

Task  2:  Analysis 

The  purpose  of  the  analysis  effort  is  to  optimize  the  algorithm  and  to 
predict  its  performance  in  term  of  probabilities  of  detection  and  of  false 
alarm.  Optimization  of  the  algorithm  will  include  tradeoffs  between 
performance  and  computational  cost. 

Task  3:  Simulation 

The  purpose  of  the  simulation  work  is  to  verify  and  complement  the 
analysis.  Each  of  the  analytical  results  need  to  be  validated  by 
simulation.  Also,  various  cases  whose  performance  is  difficult  to  analyze, 
need  to  be  evaluated  by  simulation.  An  example  of  such  a  case  is  that  of 
targets  which  exit  or  enter  the  relevant  ocean  volume  during  the  current  K- 
frame  processing  interval. 
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Appendix  A 


In  this  appendix,  the  result  of  the  CtrlC  "BROWSE"  command  is  listed. 
The  BROWSE  command  lists  all  major  functions  in  the  CtrlC  environment 
organized  by  topic. 


A. 1.0  CTRL-C  INFORMATION  AON  ASSISTANCE 


browse  -  displays  this  CTRL-C  command  synopsis  list 
help  -  help  facility 

intro  -  a  quick  introduction  for  the  new  user 
logo  -  graphs  CTRL-C  logo 

-  news  on  updates  and  new  additions 

-  lists  current  User-Defined  Functions  and  libraries 

-  lists  current  variables 

-  gives  succinct  answers  to  any  questions 


news 

what 

who 

why 


A. 2.0  STARTING,  INTERRUPTING,  AND  ENDING  THE  CTRL-C  SESSION 


exit  -  end  CTRL-C  session  and  save  current  workspace 

quit  -  end  CTRL-C  session  and  do  not  save  current  workspace 

resume  -  resume  CTRL-C  session  with  workspace  from  last  session 
S  -  interrupts  the  session  and  allows  a  VMS  command  to  be  typed 
tY  -  aborts  CTRL-C 

tC  -  local  CTRL-C  abort,  stops  current  command 


A, 3.0  BASICS 


A. 3.1  SETTING  THE  OPERATING  ENVIRONMENT 


char  -  character  set  replacement 
chop  -  truncate  arithmetic 

clear  -  erase  variables,  functions,  or  libraries  from  work  space 

disp  -  displays  variables  containing  text 

hard  -  specifies  the  hardcopy  device  type  for  graphics 

lines  -  sets  maximum  number  of  lines  of  a  matrix  that  are  displayed 
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|  long  -  sets  long  numerical  output  display  format 

!  page  -  clears  alphanumeric  CRT  screen 

semi  -  cause  or  suppress  printing  of  operation  results 
short  -  sets  short  numerical  output  display  format 
term  -  specifies  the  terminal  type  for  graphics 

f 

j  A. 3.2  SPECIAL  SYMBOLS 

| 

)  [  -  used  in  forming  vectors  and  matrices  (also  used  for  macros) 

]  -  see  [ 

<  -  means  "get  input  from  file"  or  less  than 

,  >  -  means  "send  output  to  file"  or  greater  than 

(  -  indicates  precedence  in  arithmetic  expressions  in  the  usual  way 

)  -  see  ( 

=  -  used  in  assignment  statements 

.  -  Decimal  point.  Element-by-element  operations.  Kronecker  products. 

..  -  continuation  of  a  statement  onto  the  next  line 
,  -  separates  matrix  subscripts  and  function  arguments 

;  -  Ends  rows.  Suppresses  printing  of  commands. 

/  -  matrix  right  division  computed  by  Gaussian  elimination 
//  -  comment  statement,  function  definition 
=  -  matrix  left  division  computed  by  Gaussian  elimination 
'  -  matrix  transpose,  quote  to  delimit  character  strings 

+  -  matrix  addition 

-  -  matrix  subtraction 

*  -  matrix  multiplication 
**  -  raises  matrices  to  powers 
:  -  used  in  subscripts  and  FOR  iterations 

S  -  interrupts  the  session  and  allows  a  VMS  command  to  be  typed 

A. 3. 3  SPECIAL  PERMANENT  VARIABLES 

ans  -  variable  created  when  expressions  are  not  assigned 

eps  -  floating  point  relative  accuracy 

flop  -  count  of  floating  point  operations 

hard  -  specifies  the  hardcopy  device  type  for  graphics 
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pi  -  3.1415926536 

term  -  specifies  the  terminal  type  for  graphics 
A. 3.4  PROGRAM  FLOW  STATEMENTS 


else  -  used  with  if 

end  -  delineates  'for1,  'while1,  and  'if'  statements 

exit  -  terminates  a  'for'  or  'while'  loop;  or  ends  session 

for  -  repeat  statements  a  specific  number  of  times 

if  -  conditionally  execute  statements 

while  -  repeat  statements  while  an  expression  is  true 

A. 4.0  GRAPHICS 


awi  n 
di  sp 

erase  - 
page 
plot 
pline  - 
pmark  - 
p3d 

replot  - 
redterm- 
redhard- 
text 

title  - 
window  - 
xlabel  - 
ylabel  - 


split  screen  or  window  on  terminal  in  alphanumeric  mode 

compact  matrix  print  as  +,  -,  and  blank 

erases  entire  plotting  surface 

clears  alphanumeric  CRT  screen 

X-Y  point,  bar  and  line  plots 

GKS  polyline  primitive 

GKS  polymarker  primitive 

3-dimensional  plotting  of  a  matrix 

repeat  plot  to  a  hardcopy  file 

redirect  terminal  graphics  to  a  file 

redirect  hardcopy  graphics  to  a  file 

GKS  text  placement  primitive 

plot  titles  for  PLOT  and  P3D 

specifies  the  window  for  subsequent  graphics  commands 
X-axis  labels  for  PLOT 
Y-axis  labels  for  PLOT 


A. 5.0  PROCEDURES  AND  USER-DEFINED  FUNCTIONS 


ascii  -  converts  variables  from  decimal  to  ASCII 
base  -  converts  numbers  to  different  bases 
echo  -  echoing  of  the  commands 
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emode  -  stops  CTRL-C  if  an  error  occurs 

do  -  runs  a  CTRL-C  procedure  located  in  a  file 

disp  -  displays  variables  as  character  strings 

edf  -  edit  a  function  using  the  EDT  editor 

deff  -  define  a  procedure  to  be  a  function 

glob  -  defines  global  variables  within  functions 

input  -  prompts  user  for  input 

str  -  Number  to  string  conversion 

lib  -  function  library  maintainer 

pause  -  pause  until  <cr>  entered 

return-  returns  from  a  user  function 

A. 6.0  DISK  FILES 

diary  -  saves  a  “diary"  of  the  session  in  a  disk  file 

key  -  saves  all  keyboard  input  in  a  file 

load  -  loads  variables  from  a  file 

print  -  prints  variables  on  a  file  with  132  columns 

save  -  saves  variables  on  a  file 

A. 7.0  MATRIX  ANALYSIS 

A. 7.1  E1EHENT  BY  ELEMENT  MATRIX  FUNCTIONS 

abs  -  absolute  value  or  magnitude 

conj  -  complex  conjugate 

imag  -  imaginary  part 

max  -  largest  value 

min  -  smallest  value 

prod  -  product  of  elements 

real  -  real  part 

round  -  round  to  nearest  integer 

-  sum  of  the  elements  of  a  matrix 


sum 


A. 7. 2  SQUARE  MATRIX  ELEMENTARY  FUNCTIONS 

atan  -  arctangent 

cos  -  cosine 

exp  -  matrix  exponential 

log  -  natural  logarithm 

inv  -  inverse  of  a  matrix 

sin  -  sine 

sqrt  -  square  root 

A. 7. 3  MATRIX  PROPERTIES 

cond  -  condition  number  in  2-norm 
det  -  determinant 

norm  -  singular  values,  1-norm,  infinity  norm,  and  F-norm 

rank  -  rank  of  a  matrix 

rat  -  remove  roundoff  error 

rcond  -  estimate  of  the  condition  of  a  matrix 

A. 7. 4  MATRIX  DECOMPOSITIONS  AND  FACTORIZATIONS 

chol  -  Cholesky  factorization 

eig  -  eigenvalues  and  eigenvectors 

geig  -  generalized  eigenvalues  using  the  QZ  algorithm 

hess  -  Hessenberg  form 

lu  -  factors  from  Gaussian  elimination 

orth  -  orthogonal izati on 

qr  -  orthogonal -tri angul ar  decomposition 

schur  -  Schur  decomposition 

rref  -  reduced  row  echelon  form  of  a  rectangular  matrix 
svd  -  singular  value  decomposition 

A. 7. 5  MATRIX  DIFFERENTIAL  EQUATION  SOLUTIONS 

are  -  algebraic  Riccati  equation  solution 

dare  -  discrete  algebraic  Riccati  equation  solution 


lyap  -  Lyapunov  equation  solution 

dlyap  -  discrete  Lyapunov  equation  solution 

A. 7. 6  OTHER  OPERATIONS 

conv  -  convolution  and  polynomial  multiplication 

deconv-  deconvolution  and  polynomial  division 

pinv  -  pseudoinverse 

poly  -  characteristic  polynomial 

roots  -  find  polynomial  roots 

kron  -  Kronecker  tensor  product 

size  -  row  and  column  dimensions  of  a  matrix 

sort  -  sorting 

less  -  less  than  or  greater  than 

A. 7.7  MATRIX  BUILDING  FUNCTIONS 

diag  -  puts  in/  takes  out  vectors  on  matrix  diagonals 

eye  -  generates  identity  matrices 

hi lb  -  generates  an  inverse  Hilbert  matrix 

magic  -  generates  a  magic  square 

ones  -  generates  a  matrix  of  ones 

rand  -  generates  random  numbers  and  matrices 

tril  -  lower  triangular  part 

triu  -  upper  triangular  part 

A. 8.0  CONTROL  DESIGN  AND  ANALYSIS 

A. 8.1  CONTINUOUS-TIME  SYSTEMS  DESIGN 

s  =  are(a,b,q,r)  algebraic  Riccati  equation  solution 

[k,p]  =  1qe(a,b,c,q,r)  optimal  LQG  estimator  design 

[k,s]  =  lqr(a,b,q,r)  optimal  LQG  regulator  design 

[k,s]  =  lqry(a,b,c,d,q,r,n)  LQG  regulator  design  with  output  weighting 

[k,s]  =  impl (a,b,c,f ,q,r)  implicit  model  following  gain  calculation 
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V 

k 

= 

place(a,b,p) 

regulator  and  estimator  pole  placement 

(SISO) 

1 

r 

rloc(n,d,k) 

root  loci  Analysis  Time  domain 

y 

s 

impul (a,b,c,iu,t) 

impulse  response 

V\ 

V. 

y 

= 

step(a,b,c,d,iu,t) 

step  response 

y 

= 

ramp(a,b,c,d,iu,t) 

ramp  response 

f- 

y 

= 

simu(a,b,c,d,u,t) 

simulation  with  arbitrary  inputs  Frequency 

A' 

domai n 

■v 

[m,p] 

= 

bode{a,b,c,d,iu,w) 

Bode  and  Nichols  frequency  response 

f  * 

O, 

i] 

= 

nyqu(a,b,c,d,iu,w) 

Nyquist  frequency  response 

9 

= 

freq(a,b,c,d,s) 

complex  frequency  response 

y 

sv 

= 

sigma(a,b,c,d,w) 

singular  value  frequency  response  (sigma 

plot) 

w 

= 

logspace(dl,d2) 

frequency  vector  generation  Conversions 

c2d  -  conti nous  to  discrete  state  space 
c2dt  -  continuous  to  discrete  state  space  with  pure  delay 
psit  -  intermediate  discrete  conversion  function 
ss2tf  -  state  space  to  Laplace  tranfer  function 
tf2ss  -  Laplace  transfer  function  to  state-space 

A. 8. 2  DISCRETE-TIME  SYSTEMS  DESIGN 


s 

U,pl 

[k,s] 

x 

r 

y 

[m,p] 

9 

w 

ss2tf 

tf2ss 


dare(a,b ,q,r)  discrete  algebraic  Riccati  equation 

dlqe(a,b,c,q,r)  optimal  LQG  estimator  design 

dlqr(a,b,q,r)  optimal  LQG  regulator  design 

dlyap(f,q)  discrete  Lyapunov  equation  solution 

rloc(n,d,k)  root  loci  Analysis  Time  domain 

dsimu(phi ,gam,c,d,u)  simulation  with  arbitrary  inputs 

Frequency  domain 

dbode(a,b,c,d,iu,w)  Bode  and  Nichols  frequency  response 

freq(a,b,c,d,s)  complex  frequency  response 

logspace( dl,pi )  frequency  vector  generation  Conversions 


dbode(a,b,c,d,iu,w)  Bode 
freq(a,b,c,d,s)  comp 
logspace( dl.pi)  freqi 
state  space  to  Z-transform 
Z-transform  to  state-space 
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A. 8. 3  OTHER  BUILDING  STATE-SPACE  SYSTSEMS 


series  -  series  connection 
para!  -  parallel  connection 
interc  -  block  diagram  system  interconnection 
System  Matrix  properties 
ctrb  -  controllability  staircase  form 
obsv  -  observability  staircase  form 
minreal  -  minimum  realization 
balreal  -  balanced  realization 
gram  -  controllability  and  observability  gramians 
stair  -  staircase  algorithm 
tzero  -  transmission  zeros 

A. 9.0  DIGITAL  SIGNAL  PROCESSING 

conv  -  convolution 

deconv  -  deconvolution 

fft  -  1  and  2  dimensional  FFT 

ifft  -  l  and  2  dimensional  Inverse  FFT 

tdlf  -  tapped-del ay-line  filter 

maxlike  -  maximum  likelihood  identification 

rml  -  recursive  maximum  likelihood  for  ARMAX  model 

Typing  'BROWSE  >filename'  will  send  this  listing  to  a  file  where  you  may 
print  it  on  your  local  hardcopy  device. 


Appendix  B 


In  this  appendix,  the  currently  existing  version  of  CtrlC  is  used  in  a 
simple  filter  design  exercise: 

//  Example:  Design  of  an  FIR  bandpass  filter  using  windowing 

//  Step  1:  Specify  desired  frequency  response  over  the  first 
//  half  of  a  512  point  sequence: 

Hfd  =  [0*ones( 1,25)  ones(l,20)  0*ones(l,212)];  //  [000011110000]  (Bandpass) 

//  Step  2:  Reflect  about  point  1  (mod  512)  for  the  remainder  of  the 
//  sequence.  The  frequency  specification  must  be  symmetric 

//  in  order  to  have  real  filter  coefficients. 

Hf  =  [Hfd  Hf d( 256 :-l:2)3;  //  Oesired  magnitude  frequency  response 

Ht  =  ifft(Hf);  //  Step  3:  Create  corresponding  impulse  response 

//  Step  4:  Window  the  impulse  response  for  shorter  filter  length 

NN  =  50;n  =  l:nn; 

w  =  0.5*(ones(n)  -  cos(2*pi*n/(NN-l ) ) ) ;  //  Use  a  Hanning  window 

b  =  [ ht( 488: 512 )  ht(l:25)].*w; 

//  Step  5:  Find  actual  frequency  response 

ht  =  [b  0*ones( 1 ,512-NN ) ] ; 

Hf2  =  fft(ht);  //  True  complex  frequency  response 

Overlay(Ob(Hf ) ,Db(Hf2) );  //  Plot  overlay  of  desired  and  obtained  F.R. 
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digital  signal  processing.  His  most  recent  work  has  been  in  the  areas  of 
adaptive  filtering,  system  identification,  digital  filter  design,  delay 
estimation,  spectrum  analysis,  and  signal-processing  system  development. 

Background 

Dr.  Smith  entered  Rice  University  on  a  Moody  Foundation  scholarship,  and, 
while  at  Rice,  received  the  Brown  Engineering  Merit  Award  in  two  successive 
years.  In  his  Junior  year  he  was  elected  to  Tau  Beta  Pi.  In  1975  he  finished 
his  Bachelor's  degree  (Magna  Cum  Laude)  in  electrical  engineering, 
specializing  in  control,  communications,  and  circuits.  Julius  was  fortunate 
to  have  taken  every  course  taught  by  Prof.'s  C.  S.  Burrus  and  J.  B.  Pearson. 
This  provided  him  with  a  strong  foundation  for  further  research  in  digital 
signal  processing. 

During  his  last  year  at  Rice,  Julius  worked  for  the  IBM  Houston 
Scientific  Center  codesigning  and  implementing  a  user  interface  to  a 
comprehensive  set  of  tools  for  power  systems  analysis.  The  software  package 
was  published  under  the  acronym  MAPS  for  Modular  Analysis  of  Power  Systems. 

Immediately  after  graduation  from  Rice  in  1975,  Dr.  Smith  began  working 
as  a  signal  analyst  in  training  at  the  Electromagnetic  Systems  Laboratories 
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(ESI)  in  Sunnyvale  California. 
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Or.  Smith's  first  project  at  ESL  was  to  design  and  simulate  the  digital 
AM/FM  subsystem  of  an  automatic  signal  recognition  system.  He  was  responsible 
for  high-level  hardware  specifications  and  system  performance  evaluation.  In 
connected  tasks,  he  developed  new  signal  detection  techniques  to  exploit 
advantages  inherent  in  his  AM/FM  receiver  design.  During  this  time  Dr.  Smith 
acquired  extensive  experience  in  applied  digital  signal  processing. 


In  1977,  Dr.  Smith  was  awarded  a  Hertz  Foundation  Graduate  Fellowship, 
and  in  the  Fall  of  that  year  he  became  a  full-time  graduate  student  at 
Stanford  in  the  Electrical  Engineering  Department.  The  first  year  was  devoted 
to  completing  course  work  for  the  Master's  degree,  and  he  received  his  M.S.  in 
1978  with  a  near  perfect  GPA. 

In  the  summer  of  1978  Dr.  Smith  returned  to  ESL  to  implement  the  dynamic¬ 
programming  component  of  an  automatic  language  discriminator.  The 
di scriminator  produced  a  forced  decision  among  five  different  languages  based 
on  the  likelihood  of  the  estimated  phoneme  transition  probabilities. 

Also  in  the  summer  of  1978,  Dr.  Smith  performed  consulting  work  for  Total 
Technology  Inc.  His  task  was  to  model  optical  transfer  functions  of  phosphor 
screens  used  in  image-intensifier  tubes.  This  was  Dr.  Smith's  first 
involvement  with  the  field  of  system  identification  which  became  a  major 
emphasis  of  his  thesis  work. 

In  the  summer  of  1979,  Dr.  Smith  worked  in  the  Acoustics  Research 
Department  of  Bell  Laboratories  in  Murray  Hill  Mew  Jersey.  With  Jont  B. 

Allen,  he  designed,  implemented,  evaluated,  and  published  a  new  technique  for 
adaptive  delta  modulation  (ADM).  The  new  ADM  coder  operated  with  a  variable 
sampling  rate  which  was  slaved  to  the  measured  input  signal  bandwidth. 

In  the  summer  of  1980,  Dr.  Smith  returned  to  ESL  to  carry  out  research 
applying  linear  prediction  techniques  to  FSK  demodulation.  Two  novel 
algorithms  were  developed,  evaluated,  and  documented  in  an  internal  final 
report. 
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Also  in  the  summer  of  1980,  Dr.  Smith  performed  consulting  services  for 
Narendra  Gupta  of  Integrated  Systems  Inc.  His  task  was  to  implement  a 
maximum-likelihood  system  identification  program  which  was  based  on  a  state- 
space  model  explicitly  parametrized  in  terms  of  second-order  mode  frequencies,  y 

bandwidths,  gains,  and  initial  conditions.  •-< 


In  the  summer  of  1980,  Dr.  Smith  took  a  brief  visiting  position  at  the 
Computer  Audio  Research  Laboratory,  University  of  California,  San  Diego. 

There  he  installed  the  IEEE  Programs  for  Digital  Signal  Processing  on  a  Vax 
11/780  mainframe  computer  and  produced  user-interface  software  in  the  C 
programming  language  which  made  extensive  use  of  the  Unix  pioes  facility.  The 
result  was  an  elegant  and  powerful  set  of  signal  processing  tools  which 
continue  to  be  widely  distributed  by  UCSD  in  the  public  domain. 

By  1980  Dr.  Smith  had  completed  his  Ph.D.  course  work  at  Stanford,  and  he 
devoted  the  next  two  years  to  his  thesis  research.  The  first  area  of  intense 
study  was  system  identification.  During  this  time  Prof.  Lennart  Ljung  was  a 
visiting  professor  at  Stanford,  and  Dr.  Smith  benefited  greatly  from  pursuing 
problems  under  Prof.  Ljung' s  direction.  In  addition  to  developing  techniques 
for  system  identification,  Dr.  Smith  developed  several  contributions  to 
digital  filter  design  (a  closely  related  area).  His  work  in  both  areas 
concentrated  on  effective  manipulation  of  the  error  criterion  minimized  by 
modeling  techniques.  Some  of  the  most  useful  improvements  were  obtained 
through  the  use  of  conformal  mapping  techniques  and  other  preprocessing 
steps.  He  proved  that  there  is  no  upper  bound  on  the  number  of  false  local 
minima  when  miiimizing  output  error  with  respect  to  even  a  single  filter 
pole.  His  thesis  contains  what  is  perhaps  the  first  use  of  Hankel  norm  theory 
to  obtain  a  globally  convergent  algorithm  for  optimally  approximating  a 
complex  desired  frequency  response  with  a  recursive  digital  filter.  The  Ph.D. 
degree  was  awarded  in  June  1983.  Later  that  year  he  was  elected  to  the  Sigma 
Xi  society. 

In  the  Fall  of  1982,  Dr.  Smith  joined  the  Adaptive  Systems  Department  of 
the  Advanced  Technology  Division  of  Systems  Control  Technology  (SCT)  in  Palo 
Alto  California.  He  worked  half  time  until  December  1982  whereupon  he  became 
a  full  time  SCT  employee.  He  continued  full  time  until  October  of  1984  at 


v. 
*  » 

L3 


50 


which  time  he  reduced  his  commitment  to  60  percent  in  order  to  make  time  to 
teach  a  two-year  course  sequence  in  digital  signal  processing  at  Stanford. 

Dr.  Smith  has  worked  exclusively  on  projects  directed  by  Dr.  Benjamin 
Friedlander  at  SCT  for  the  past  three  years. 

Tasks  at  SCT 

The  overall  theme  of  recent  research  in  the  Apdaptive  Systems  Department 
at  SCT  has  been  to  develop  modern  algorithms  which  address  problems 
fundamental  to  target  localization  based  on  multiple  hydrophone  recordings. 

The  principal  new  technologies  brought  to  bear  on  the  localization  objective 
have  been  drawn  from  the  most  recent  advances  in  system  identification  and 
parametric  spectral  estimation.  Complementing  this  line  of  research,  Dr.  Smith 
has  developed  more  traditional  sonar  signal  processing. facilities  for 
nonparametric  spectral  estimation  based  on  the  FFT. 

The  first  project  assigned  to  Dr.  Smith  in  1983  was  to  develop  an 
adaptive  time-delay  estimation  algorithm  using  fundamental  formulations  from 
recursive  system  identification.  This  effort  resulted  in  a  highly  efficient 
and  accurate  method  with  a  very  rapid  tracking  capability  in  the  time-varying 
case.  The  heart  of  the  method  was  a  novel  intersample  interpolator  developed 
by  Dr.  Smith  for  his  thesis  applications. 

The  adaptive  delay  estimation  technique  was  next  extended  in  two 
straightforward  ways  to  the  problem  of  adaptively  tracking  multipath  delay. 

The  techniques  were  compared  in  a  precise  way  to  the  maximum  likelihood 
estimator  for  the  multipath  estimation  problem.  Theoretical  characterizations 
of  the  asymptotic  properties  of  the  estimates  were  carried  out,  including 
global  convergence  studies,  derivation  of  asymptotic  parameter  variance  (with 
comparisons  to  the  Cramer-Rao  lower  bound),  analysis  of  multipath 
detection/false-alarm  probabilities,  and  the  effects  of  multipath  on  target 
local i zation. 

A  third  task  was  to  develop  an  adaptive  notch  filter  algorithm,  again 
applying  the  latest  results  in  recursive  identification  as  well  as  some  new 
theorems  on  rational  transfer  functions  developed  by  Dr.  Smith  in  the  course 


of  his  thesis  work. 


A  task  carried  out  in  parallel  with  the  above  identification  projects  was 
the  ongoing  evolution  of  the  Modified  Yule-Walker  (MYW)  method  for  ARMA 
spectral  estimation.  More  recent  variations  of  the  MYW  approach  have  been 
termed  Optimal  Instrumental  Variable  methods.  The  MYVI  methods  can  be  viewed 
as  variations  of  Prony's  method  for  pole-zero  digital  filter  design  using 
estimates  of  the  correlation  function  in  place  of  the  filter  impulse  response. 

The  numerous  MYW  enhancements  were  aimed  at  achieving  reduced  variance 
and  bias  in  the  parameter  estimates  for  a  given  number  of  data  observations. 
Dr.  Smith  performed  all  simulations  and  software  development  for  these  studies 
and  contributed  major  implementation  decisions.  Very  significant  improvements 
over  the  original  MYW  technique  were  realized.  To  aid  in  the  dissemination  of 
these  superior  techniques.  Dr.  Smith  wrote  highly  portable  Fortran  versions  of 
the  MYW  methods. 

A  later  task  was  to  implement  a  technique  for  multiple  target 
localization  from  narrowband  Doppler  measurements.  The  method  consists  of 
finding  intersections  between  hypersurfaces  in  a  five-dimensional  space.  Each 
target  appears  (in  the  noiseless  case)  as  a  point  of  intersection  of  all  the 
surfaces.  Dr.  Smith  developed  novel  means  of  analyzing  the  behavior  of  this 
algorithm  in  the  presence  of  noisy  measurements.  Dr.  Smith  also  developed  a 
clustering  algorithm  for  the  technique,  including  a  method  for  determining  the 
number  of  targets,  associating  spectral  lines  to  targets,  and  assigning  points 
of  intersection  to  targets.  As  in  all  projects  on  which  he  worked.  Dr.  Smith 
was  responsible  for  all  software  development  and  computer  simulations. 

In  another  parallel  series  of  projects,  Dr.  Smith  developed  recursive 
identification  counterparts  to  the  Constant  Modulus  Algorithm  (CMA)  which  is  a 
novel  method  of  adaptive  channel  equalization  pioneered  by  J.  R.  Treichler. 

Dr.  Smith's  enhancements  included  extension  to  a  nonrecursive  channel  model 
(ideal  for  removing  multipath  distortion),  a  real-signal  version  with 
convergence  properties  as  good  as  the  complex-signal  version,  a  quadratic 
Newton-method  to  replace  fixed-step-size  gradient  descent,  an  interference 
rejection  technique,  and  more.  In  addition,  Dr.  Smith  first  proved  global 


convergence  of  the  CMA  from  an  arbitrary  starting  point  for  the  model -complete 
case.  This  was  perhaps  the  first  global  convergence  proof  in  the  system 
identification  literature  which  involved  a  prediction  error  quadratic  rather 
than  linear  in  the  model  parameters.  The  notion  of  "persistently  exciting" 
was  generalized  to  the  nonsingularity  of  a  rank  four  "covariance  tensor." 

To  obtain  useful  benchmarks  of  comparison  between  the  ARMA  spectral 
estimation  methods  (MYW)  and  the  perceived  standard  practice,  Dr.  Smith 
developed  an  implementation  of  the  periodogram  method  for  spectrum  estimation 
(PMPSE)  using  the  FPS  array  processor.  The  PMPSE  essentially  averages  the 
squared  magnitude  of  successive  overlapped,  Hamming-windowed,  short-time 
FFT's.  In  an  internal  report  by  Dr.  Smith  containing  the  MYW  and  PMPSE 
programs,  a  review  of  the  theoretical  asymptotic  performance  of  both 
approaches  was  included.  Dr.  Smith  also  developed  a  wide  variety  of  spectral 
displays,  including  ordinary  line  plots,  "waterfall  diagrams,"  encoded  gray- 
level  laser-printer  plots,  and  he  investigated  the  benefits  of  using  edge- 
enhancement  algorithms  already  available  at  SCT. 

To  allow  benchmark  comparisons  between  the  new  techniques  developed  by 
Dr.'s  Friedlander  and  Smith  and  previously  existing  localization  systems.  Dr. 
Smith  is  currently  finishing  up  a  replication  of  what  is  considered  standard 
low-level  sonar  signal  processing  practice.  Estimates  of  intersensor  time 
delay,  multipath  time  delay.,  intersensor  differential  Doppler,  and  absolute 
Doppler  are  all  derived  in  one  way  or  another  from  the  FFT-based  spectrum 
estimate  (using  the  PMPSE).  Both  narrowband  and  broadband  Doppler  estimation 
techniques  have  been  developed.  Excellent  differential  and  absolute  Doppler 
estimates  have  been  produced  from  live  sonar  recordings. 
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YAIR  BARNIV,  Senior  Engineer 


Dr.  Barniv  received  his  B.S.  degree  in  Electrical  Engineering  from 
Israel-Institute-of-Technology  in  1963.  He  received  his  M.S.  and  Ph.D. 
degrees  in  communications,  stochastic  processes,  and  signal/image  processing 
from  Carnegi e-Mellon  University  (CMU)  in  1978  and  1981,  respectively.  His 
Ph.D.  thesis  work  was  on  "Multi-Sensor  Image  Registration". 


Dr.  Barniv  served  as  a  commissioned  officer  in  the  Israeli  Air-Force  from 
June  1963  until  May  1968,  working  in  maintenance  and  development  of  airborne 
and  ground  radars,  and  fire  control  systems.  He  was  also  responsible  for  air 
traffic  and  automatic  landing  instrumentation  (GCA). 


From  May  1968  to  September  1977,  Dr.  Barniv  was  with  Israel  Aircraft 
Industries  (IAI).  Until  January  1970,  he  served  as  an  electrical  and  control 
engineer  designing  a  ground-based  navigation  system.  From  January  1970  to 
April  1971,  he  worked  as  a  control  engineer  and  designed  the  control  loops  for 
a  missile,  including  the  autopilot,  the  guidance,  and  the  homing  laws  and 
logic.  He  also  conducted  air-tunnel  measurements  and  the  necessary 
aerodynamic  data  extraction. 


Between  April  1971  and  January  1974,  Dr.  Barniv  was  with  a  Systems 
Engineering  Group,  whose  tasks  included  performance  estimation,  feasibility, 
and  cost-effectiveness  studies  in  fields  such  as  electro-acoustics,  electro¬ 
optics,  radar,  and  navigation  systems. 


From  1974  to  May  1975  he  served  as  a  senior  systems  engineer,  heading  a 
group  tasked  with  preparing  overall  specifications  for  a  missive  design 
project,  and  coordinating  the  design  and  implementation  activities  and  final 
integration.  He  developed  the  simulation  programs  and  used  them  to  optimize 
the  design  parameters  and  define  the  limits  of  performance.  He  was  also 
responsible  for  the  planning  and  management  of  the  flight  test  program  and 
evaluation. 


From  May  1975  to  September  1977,  he  headed  the  Systems  Engineering  Group 
working  on  the  control  issues  of  the  Mass  Transit  System  for  the  City  of 
Tel -Aviv.  His  task  was  to  study  existing  systems  from  the  point  of  view  of 
automatic  train  control,  operation,  and  supervision;  this  study  included  the 
issues  of  communications  and  control  and  man-machine  interfaces. 


From  September  1977  Dr.  Barniv  worked  on  his  Ph.D.  thesis  at  Carnegie- 
Mellon  University,  and  in  May  1979  joined  SCI  as  a  research  engineer.  He  is 
currently  working  in  the  areas  of  syntheti c-aperture-radar,  flight  path 
optimization,  closed-loop  ECM  techniques  against  various  types  of  radars, 
missile's  mi ss-oi stance  optimization,  and  dynamic  programming  application  to 
mosaic  sensor  tracking. 
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Dr.  Barnlv's  publications  Include: 

1.  "A  Ground-based  Navigation  System,"  Internal  IAI,  January  1970, 
(CONFIDENTIAL). 

2.  "A  Compensation  Method  for  the  Random  Effect  in  a  Proportional  Navigation 
System,"  internal  IAI  February  1971,  (CONFIDENTIAL). 

3.  "A  Real-time  G-Bias  Compensation  Method,"  internal  IAI,  August  1970, 
(CONFIDENTIAL). 

4.  "Location  of  P-  is  Through  the  Use  of  Acoustical  Detection,"  internal  IAI, 
October  1971,  (CONFIDENTIAL). 

5.  "Time  of  Arrival  (TOA)  Navigation,"  internal  IAI,  May  1972, 
(CONFIDENTIAL). 

6.  "The  Effects  of  the  Electrostatic  Charging  of  an  Aircraft  on  the  On-board 
Instrumentation,"  internal  IAI,  March  1973,  (CONFIDENTIAL). 

7.  "Automatic  Air  Raid  of  Naval  Targets,"  internal  IAI,  October  1973, 
(CONFIDENTIAL). 

8.  "A  Simple  Algorithm  for  Fire-Zone  Determination  for  a  Ground- to-Ground 
Missile,"  internal  IAI,  November  1974  (CONFIDENTIAL). 

9.  "A  Centralized  Supervision  System  for  a  Mass-Transit  System,"  internal 
IAI,  December  1976. 

10.  "Intensity  Disparity  Compensation  in  Optical  Pattern  Recognition,"  M.S. 
Thesis,  Carnegi e-Mel Ion  University,  May  1978. 

11.  "Correlation  of  Images  with  Random  Contrast  Reversals,"  (with  H, 

Mostafavi  and  David  Casasent),  SPIE  Vol.  238,  Image  Processing  for 
Missile  Guidance,  p.  156,  1980. 
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